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ABSTRACT 

Many extrasolar planets orbit sufficiently close to their host stars that sig- 
nificant tidal interactions can be expected, resulting in an evolution of the spin 
and orbital properties of the planets. The accompanying dissipation of energy 
can also be an important source of heat, leading to the infiation of short-period 
planets and even mass loss through Roche-lobe overfiow. Tides may therefore 
play an important role in determining the observed distributions of mass, orbital 
period, and eccentricity of the extrasolar planets. In addition, tidal interactions 
between gaseous giant planets in the solar system and their moons are thought to 
be responsible for the orbital migration of the satellites, leading to their capture 
into resonant configurations. 

Traditionally, the efficiency of tidal dissipation is simply parametrized by a 
quality factor Q, which depends, in principle, in an unknown way on the frequency 
and amplitude of the tidal forcing. In this paper, we treat the underlying fiuid 
dynamical problem with the aim of determining the efficiency of tidal dissipation 
in gaseous giant planets such as Jupiter, Saturn, or the short-period extrasolar 
planets. Efficient convection enforces a nearly adiabatic stratification in these 
bodies, which may or may not contain rocky cores. With some modifications, 
our approach can also be applied to fully convective low-mass stars. 

In cases of interest, the tidal forcing frequencies are typically comparable to 
the spin frequency of the planet but are small compared to its dynamical fre- 
quency. We therefore study the linearized response of a slowly and possibly 
differentially rotating planet to low-frequency tidal forcing. Convective regions 
of the planet support inertial waves, which possess a dense or continuous fre- 
quency spectrum in the absence of viscosity, while any radiative regions support 
generalized Hough waves. We formulate the relevant equations for studying the 
excitation of these disturbances and present a set of illustrative numerical calcu- 
lations of the tidal dissipation rate. 
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We argue that inertial waves provide a natural avenue for efficient tidal dissi- 
pation in most cases of interest. The resulting value of Q depends, in principle, 
in a highly erratic way on the forcing frequency, but we provide analytical and 
numerical evidence that the relevant frequency-averaged dissipation rate may be 
asymptotically independent of the viscosity in the limit of small Ekman number. 
For a smaller viscosity, the tidal disturbance has a finer spatial structure and 
individual resonances are more pronounced. 

In short-period extrasolar planets, tidal dissipation via inertial waves becomes 
somewhat less efficient once they are spun down to a synchronous state. However, 
if the stellar irradiation of the planet leads to the formation of a radiative outer 
layer that supports generalized Hough modes, the tidal dissipation rate can be 
enhanced, albeit with significant uncertainty, through the excitation and damping 
of these waves. The dissipative mechanisms that wc describe offer a promising 
explanation of the historical evolution and current state of the Galilean satellites 
as well as the observed circularization of the orbits of short-period extrasolar 
planets. 

Subject headings: hydrodynamics — planets and satellites: general — waves 

1. INTRODUCTION 

1.1. Tidal interactions in planetary systems 

The discovery of the very first extrasolar planet orbiting a main-sequence star, 51 Peg b, 
brought the surprising revelation that other planetary worlds can revolve around their host 
stars in extraordinarily tight orbits (Mayor & Queloz 1995). Around approximately one 
per cent of all the stars targeted by radial-velocity searches, similar planets are found with 
periods of a few days. These short-period planets were probably formed within a protostellar 
disk several AU away from their host stars through a sequence of physical processes outlined 
in the conventional 'core-accretion' theories of planetary formation (Lissauer 1993; Pollack 
et al. 1996). During its final accretion phase, strong interactions between a protoplanet 
and its natal disk lead to the formation of a gap near its orbit, which effectively terminates 
its growth. Angular momentum exchange with the disk causes a protoplanet formed in the 
inner region of the disk to migrate towards its host star (Goldreich & Tremaine 1980; Lin 
& Papaloizou 1986). Near the star, the protoplanct's migration may be halted either by 
its tidal interaction with a rapidly spinning host star or by entering a cavity in the disk 
associated with the stellar magnetosphere (Lin, Bodenheimer, & Richardson 1996). Planets 
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may also terminate their migration with intermediate-period orbits when their natal disks 
are rapidly depleted (Trilling, Lunine, & Benz 2002; Armitage et al. 2002; Ida & Lin 2003). 

Short-period planets may be scattered to the vicinity of (or far away from) the stel- 
lar surface owing to dynamical instabilities associated with their long-term post-formation 
gravitational interaction (Rasio & Ford 1996; WeidenschiUing Sz Marzari 1996). In this case, 
the protoplanet's orbit may initially be highly eccentric and subsequently undergo circular- 
ization owing to the tidal dissipation within its envelope (Rasio et al. 1996). Scattering by 
residual planetesimals is another promising avenue for planetary orbital migration (Murray 
et al. 1998). It has also been suggested that short-period planets may be formed in situ 
through a sequence of planetesimal migration, coagulation, and gas accretion (Papaloizou & 
Terquem 1999; Bodenheimer, Hubickyj, & Lissauer 2000; Sasselov & Lecar 2000). Finally, a 
totally different theory of protoplanetary formation relies upon gravitational instability in a 
massive gaseous protostellar disk (Kuiper 1951; Cameron 1978; Boss 1997), which produces 
a subcondensation with a composition close to that of the central star and having little, if 
any, solid core. Although the inferred presence of solid cores in Jupiter, Saturn, Uranus, and 
Neptune (Hubbard 1984) favors the orderly growth scenario in our solar system, there is as 
yet no direct evidence for or against the existence of cores in extrasolar planets. 

In the past eight years, more than 100 planets have been discovered, around approxi- 
mately ten per cent of the target stars in various planet search programs. A comprehensive 
distribution of planetary mass, period, and eccentricity has begun to emerge. These data are 
particularly useful for constraining and differentiating some scenarios for the formation and 
evolution of short-period planets. For example, there appears to be a minimum cut-off pe- 
riod at about three days. If the short periods of some planets were attained through orbital 
migration, their observed period distribution could be attributable to the stopping mecha- 
nism involved. It is also possible that the present period distribution is established through 
post-formation evolution, e.g., planets with extremely short periods may have perished sub- 
sequent to their formation. A particularly interesting characteristic of the planet-bearing 
stars is that they tend to be metal-rich with respect to the Sun and the F-G field-star 
average in the solar neighborhood (e.g., Gonzalez et al. 2001; Santos, Israelian, & Mayor 
2001). This correlation may result from planets having being consumed following either 
migration through the protoplanetary disk (Lin 1997; Laughlin & Adams 1997; Sandquist 
et al. 1998) or gravitational interaction with other planets (Rasio & Ford 1996). It may also 
be interpreted as evidence that an enhanced metallicity in the planet-forming disk, probably 
resulting from a metal-rich parent cloud, is especially conducive to planet formation (Pollack 
et al. 1996; Fischer & Valenti 2003; Ida & Lin 2003). 

Another important set of observational data is the planets' eccentricities. While the 
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eccentricities of planets with periods P > 21 days are uniformly distributed up to about 
0.7, all planets with P < 7 days have nearly circular orbits. This observed dichotomy in 
the eccentricity-period distribution has been attributed to the circularization of the orbits 
of short-period planets induced by their internal tidal dissipation (Rasio et al. 1996; Marcy 
et al. 1997). Such a scenario would be necessary to account for the small eccentricities of 
the short-period planets if they were scattered into, or strongly perturbed in, the vicinity of 
their host stars by their planetary siblings. If, instead, the short-period planets acquired their 
small semi-major axes through tidal interaction with their natal disks (Lin et al. 1996), such 
a process may be able to damp the eccentricities of planets with masses Mp < 10 Mj (where 
Mj is the mass of Jupiter) (Goldreich & Tremaine 1980; Artymowicz 1992), although this 
hypothesis remains uncertain (Papaloizou, Nelson, & Masset 2001; Goldreich & Sari 2003). 
However, some of the potential stopping mechanisms, such as the planets' tidal interaction 
with rapidly spinning host stars (Lin et al. 1996), can also excite their orbital eccentricities 
(Dobbs-Dixon, Lin, & Mardling 2003). An efficient post-formation eccentricity damping 
mechanism may still be needed to account for the small eccentricities of the short-period 
planets. 

Although no information is available on the rotation of extrasolar planets, tidal inter- 
action with the host star is expected to bring the rotation of a short-period planet into 
synchronism with its orbit and to eliminate any obliquity of its rotation axis. Since the 
spin of a planet accounts for much less angular momentum than its orbit, synchronization 
proceeds more rapidly than circularization. This assumption neglects the possible effects of 
thermal tides, which may tend to drive the planet away from synchronism (Thomson 1882; 
Gold & Soter 1969). 

Tidal interaction can also modify the radii and internal structures of short-period plan- 
ets. In principle, the actual sizes of extrasolar planets may be used, at least as a partial test, 
to infer whether they formed through core accretion or through gravitational instability of 

massive gaseous protostellar disks. In the fully convective envelopes of gaseous giant planets, 
it is difficult for the refractory elements to condense into droplets and become differentiated 
from the volatile elements (cf. Guillot et al. 2003). Therefore planets formed through or- 
derly core accretion are more likely to have solid cores and to be relatively compact, whereas 
those formed through gravitational instability are likely to retain a nearly uniform solar 
composition throughout their interiors and to be more extended. 

The detection of a transiting planet around the star HD 209458 provides an opportunity 
for us to measure its radius directly and thereby to constrain its present internal structure. 
The observed radius, 1.35 of this 0.63 Mj planet (Brown et al. 2001) is larger than that 
expected for a coreless planet with a similar mass and age. A planet with a core would be still 
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more compact, leading to an even larger discrepancy with the observational measurement. 
For this short-period planet, the presence of a small residual orbital eccentricity or non- 
synchronous rotation could lead to internal tidal dissipation with a heating rate comparable 
to or larger than that released by the Kelvin-Helmholtz contraction. Provided that the 
dissipation of the host star's tidal disturbance occurs well below the planet's photosphere, 
it increases the planet's internal temperature and equilibrium radius (Bodenheimer et al. 
2000). 

Planet-star interactions may also have altered the mass distribution of the short-period 
planets. Above a critical eccentricity, which is a function of the planet's semi-major axis, 
tidal dissipation of energy during the circularization process can cause a planet to inflate 
in size before its eccentricity is damped. For moderate eccentricities, the planet adjusts 
through a sequence of stable thermal equilibria in which the rate of internal tidal dissipation 
is balanced by the enhanced radiative flux associated with the planet's enlarged radius. For 
sufficiently large eccentricities, the planet swells beyond 2 Rj and its internal degeneracy is 
partially lifted. Thereafter, the thermal equilibria become unstable and the planet undergoes 
runaway inflation until its radius exceeds the Roche radius (Gu, Lin, & Bodenheimer 2003). 
The critical eccentricity of about 0.2 (for a young Jupiter-mass planet with a period less 
than three days) may be easily attained, with the result that many short-period planets may 
have migrated to the vicinity of their host stars and perished there. 

The above discussion clearly indicates the importance of planet-star tidal interactions to 
the origin and destiny of short-period planets. They may determine (1) the asymptotic semi- 
major axis of migrating protoplanets, (2) the absence of planets with periods less than three 
days, (3) the structure of those planets that survive in the vicinity of their host stars, and 
(4) the small eccentricities and relatively low masses of the short-period planets. However, 
the efficiency of these processes depends crucially on the ability of the planet to dissipate 
tidal disturbances, which is poorly understood. 

1.2. Theory of the equilibrium tide 

An early analysis of the nature of the tidal interaction between Earth and the Moon was 
introduced by Darwin (1880) based on the concept of the equihbrium tide (see Cartwright 
2000 for a historical discussion). In this model, due originally to Newton, a homogeneous 
spherical body continually adjusts to maintain a state of quasi-hydrostatic equilibrium in the 
varying gravitational potential of its orbital companion. Darwin introduced a phase lag into 
the response, the lag being proportional to the tidal forcing frequency and attributable to the 
viscosity of the body. The phase lag gives rise to a net tidal torque and dissipation of energy. 
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Subsequent formulations (Munk & MacDonald 1960; Goldreich & Soter 1966) parametrize 
the efficiency of the tidal dissipation, whatever its origin, by a specific dissipation function or 
Q- value (quality factor). Although, in principle, Q is an unknown function of the frequency 
and amplitude of the tidal forcing, in the planetary science community it is usually treated 
as a constant property of each body in the solar system, corresponding to a constant phase 
lag of arcsin(5~^ (e.g., Murray & Dermott 1999). Indeed, studies of the Earth's rotation 
provide evidence that its value is approximately constant over a wide range of frequency 
(Munk & MacDonald 1960), even though a variety of mechanisms may be responsible. The 
tidal Q- value may be predominantly determined by turbulent dissipation in the shallow seas 
(Taylor 1920; Jeffreys 1920; see Munk & MacDonald 1960 for a discussion). In the case of 
the Moon, the observational constraints on Q are much weaker, but the synchronized spin 
of the Moon suggests that the presence of an ocean is not necessary for tidal dissipation. In 
the context of binary stars, the assumption of a constant lag time is usually adopted (Hut 
1981; Eggleton, Kiseleva, & Hut 1998) as in Darwin's viscous model. 

Evidence of tidal dissipation can also be found in gaseous giant planets such as Jupiter 
and Saturn. For example, the Laplace resonance of lo, Europa, and Ganymede may be en- 
tered through the tidal interaction of the moons with Jupiter (Goldreich 1965; Peale, Cassen, 
& Reynolds 1979; Lin & Papaloizou 1979b). In order to maintain this resonant configura- 
tion despite the dissipation inside lo, tidal dissipation within Jupiter must continually induce 
angular momentum transfer from its spin to lo's orbit. The magnitude of the Q- value for 
Jupiter, inferred from lo's dissipation rate, is in the range 6 x 10^ — 2 x 10^ (Yoder & Peale 
1981). These values are similar to those inferred, for solar-type stars, from the circularization 
of close binaries as a function of their age and semi-major axis (Mathieu 1994; Terquem et al. 
1998). Using these values of Q, the orbital evolution of a planetary or a stellar companion 
of a main-sequence star can be estimated by extrapolation, under the assumption of either 
a constant lag angle (Goldreich & Soter 1966), a constant lag time (Hut 1981; Eggleton et 
al. 1998), or an intermediate approach (Mardling & Lin 2002). 

In the extended convective envelopes of gaseous giant planets and low-mass stars, turbu- 
lence can lead to dissipation of the motion that results from the continual adjustment of the 
equilibrium tide. However, the turbulent viscosity estimated from the mixing-length theory 
ought to be reduced by a frequency-dependent factor owing to the fact that the convective 
turnover timescale is usually much longer than the period of the tidal forcing. Estimat- 
ing this reduction in the efficiency of the turbulent viscosity, Zahn (1977, 1989) analyzed 
the dissipation of equilibrium tides in low-mass stars and computed an efficiency of angular 
momentum transfer for a non-rotating star that matches the observationally inferred cir- 
cularization timescale of solar- type binary stars (Mathieu 1994). However, with a different 
prescription for the efficiency reduction factor (Goldreich & Keeley 1977; Goodman & Oh 
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1997), the derived rate of angular momentum transfer falls short of that required by nearly 
two orders of magnitude (Terquem et al. 1998). Based on a similar approach, which gives 
Q ~ 5 X 10^^ from the dissipation of the equilibrium tide in Jupiter, Goldreich & Nicholson 
(1977) suggested that the tidal interaction between the Galilean satellites and Jupiter cannot 
drive any significant orbital evolution. For the short-period extrasolar planets, such a high 
Q-value would imply circularization times considerably longer than the ages of their host 
stars. 



1.3. Dynamical tides 

The problem of tidal dissipation is even more acute for high-mass stars in close bina- 
ries because they have extended radiative envelopes that may be expected to be free from 
turbulence. Yet their orbits are indeed circularized despite their short lifespans (Primini, 
Rappaport, & Joss 1977). In these systems, the tidal perturbation of the companion can 
induce the resonant excitation of low-frequency g-mode oscillations in the radiative region, 
which carry both energy and angular momentum fluxes (Cowling 1941; Zahn 1970). This 
wavelike, dynamical tide exists in addition to the equilibrium tide and provides an alternative 
avenue for tidal dissipation. The g modes are primarily excited in the radiative envelope close 
to the convective core where the Brunt-Vaisala frequency is comparable to the tidal forcing 
frequency and the wavelengths of the gravity waves are sufficiently long to couple well with 
the tidal potential (Goldreich & Nicholson 1989). Using asymptotic methods in the limit 
that the forcing frequency is small compared to the dynamical frequency of the star, Zahn 
demonstrated the existence of resonant modes that amplify the dynamical response of the 
star at particular frequencies. This local analysis was generalized to modest frequencies in 
a numerical analysis by Savonije & Papaloizou (1983, 1984). Despite the absence of convec- 
tively driven turbulence in the radiative envelope, the g-mode oscillations can be dissipated 
through radiative damping on the wavelength scale. The radiative loss in the stellar interior 
is generally weak so that the waves can propagate in an approximately adiabatic manner 
to the stellar surface where they are dissipated and the angular momentum they carry is 
deposited. This concept gives rise to the interesting suggestion that the stars become syn- 
chronized from the outside in (Goldreich & Nicholson 1989). As the dissipation layer first 
establishes a synchronized rotation, it presents a barrier to the outwardly propagating waves 
in the form of a corotation resonance at which the group velocity of the waves vanishes. 
Just below the corotation resonance, waves are stalled with large amplitude and dissipated. 
This process continues, initially inducing differential rotation but eventually leading to the 
global synchronization of the stellar spin. This theoretical model is in agreement with the 
observed spin and orbital evolution of high-mass close binary stars (Giuricin, Mardirossian, 
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& Mezzetti 1984a, 1984b). 

In the interior of a solar-type star, where a radiative core is surrounded by a convec- 
tive envelope, g-mode oscillations are also excited in the radiative region and, inasmuch 
as they extend into the outer convective region, they are dissipated by turbulent viscosity. 
For some special forcing frequencies, these g-mode oscillations can also attain a global reso- 
nance that strongly enhances the oscillation amplitude in the radiative region and the rate 
of dissipation in the convective region (Terquem et al. 1998). Nevertheless, the overall syn- 
chronization process is limited by the less efficient tidal dissipation that occurs in between 
resonant frequencies. A calibration with the observed circularization periods suggests that 
the required viscosity is nearly two orders of magnitude larger than that inferred from the 
standard mixing-length model for turbulent convection and the Goldreich & Keeley (1977) 
reduction factor. 

Although g-mode excitation is a powerful process that drives the dynamical tidal re- 
sponse in stars, it can be effective only in radiative regions. Since the envelopes of gaseous 
giant planets are mostly convective, the relevance of g-mode oscillations is less well estab- 
hshed. In a series of papers, loannou & Lindzen (1993a, 1993b, 1994) considered the dynam- 
ical response of Jupiter to the tidal perturbation of lo with a prescribed model for Jupiter's 
interior. They showed that g modes can be resonantly excited if the interior is slightly suba- 
diabatic to the extent that an appropriately tuned wave cavity is created. When these waves 
are transmitted into the atmosphere and dissipated, they can produce an effective Q- value 
comparable to that inferred from the supposed orbital evolution of lo. However, current 
models suggest that convection may extend throughout the entire envelope of Jupiter and 
is such an efficient heat-transport process that Jupiter's envelope is adiabatically stratified 
and neutrally buoyant to a very high degree (Guillot et al. 2003). In this limit, loannou 
& Lindzen's mechanism does not work and only a non-wavelike dynamical response to lo's 
tidal perturbation would remain with a much reduced amplitude and Q > 10^^. 

In contrast to Jupiter, the surface layers of short-period extrasolar planets may be 
stabilized against convection because they are intensely heated by their host stars and may 
attain a radiative state. If so, g-mode oscillations may be excited in the radiative layer just 
above its interface with the planet's convective envelope and dissipated through radiative or 
nonlinear damping (Lubow, Tout, & Livio 1997) as suggested for high- mass stars. However, 
the one-sided stellar irradiation of the surface of a synchronized planet induces a large-scale, 
shallow circulation (Burkert et al. 2003), which does not necessarily suppress convection 
and enhance the Brunt- Vaisala frequency over an extended region in the atmosphere. Thus 
the efficiency of dynamical tidal dissipation via g-mode oscillations remains an outstanding 
issue. 



-9- 



In their calculation loannou & Lindzen (1993a, 1993b) also considered the effects of 
uniform rotation. For the excitation of the dynamical response, they adopted the so-called 
'traditional approximation' in which only the radial component of the angular velocity is 
included in the computation of the Coriolis force (Eckart 1960; Chapman & Lindzen 1970; 
Unno et al. 1989). This approach was justified in their model on the basis that the dy- 
namical response occurs primarily near the atmosphere of Jupiter where the horizontal scale 
of the motion is large and the radial velocity perturbation is relatively small. It enabled 
them to separate the radial and angular variables in the governing linearized equations for 
the perturbations. The meridional structure of Jupiter's response due to Jo's perturbation 
is then determined by Laplace's tidal equation, the solutions of which are Hough modes 
(Hough 1897, 1898) instead of spherical harmonics, while the radial structure is determined 
by a set of ordinary differential equations. Using a local (WKB) approximation, they ob- 
tained a dispersion relation describing a mixture of gravity, inertial, and acoustic oscillations. 
Although inertial waves can be excited in the nearly adiabatic planetary interior, loannou 
& Lindzen (1993b) did not consider them relevant for tidal dissipation. In fact the main 
effect of rotation in their model is that tidal forcing by a single solid harmonic projects on 
to many Hough modes, each with the potential to resonate. They noted that the traditional 
approximation of neglecting the latitudinal component of the angular velocity is, in fact, not 
appropriate in the planet's interior. 

The effects of uniform rotation have been investigated in greater depth by Savonije, 
Papaloizou, & Alberts (1995) in the context of high-mass binary stars, which consist of 
extended radiative envelopes around small convective cores. These authors were the first 
to attempt a two-dimensional numerical solution of the linearized equations governing tidal 
disturbances, including the full Coriolis force. By 'switching off' terms deriving from the 
non-radial component of the angular velocity, they were able to verify that the traditional 
approximation can give a reasonably accurate measure of the tidal torque on a highly strati- 
fied, slowly rotating star. Like loannou & Lindzen (1993b), they found that the possibilities 
for resonant excitation of g modes are enhanced in a rotating star because the Coriohs force 
couples spherical harmonics of different degrees. Recently Savonije & Witte (2002) apphed 
the traditional approximation to carry out a high-resolution study of the dynamical tidal 
interaction between a uniformly rotating solar-type star and an orbital companion including 
the effects of stellar evolution. 

At the same time, Savonije et al. (1995) noted that the Coriolis force introduces a 
rich spectrum of inertial modes in a certain frequency range, which are not adequately 
represented by the traditional approximation. When the full Coriolis force is included, the 
response of the star to tidal forcing in the frequency range of the inertial modes contains large- 
amplitude, very short-wavelength components that could not be resolved on the grid. Using 
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an improved finite-difference numerical scheme including viscosity (Savonije & Papaloizou 
1997) and an asymptotic perturbation analysis (Papaloizou & Savonije 1997), they later 
showed that inertial modes can be resonantly excited in the adiabatic convective cores of 
high-mass stars and can interact with the rotationally modified g modes and torsional r modes 
in the radiative envelopes. When these waves are dissipated through radiative damping, the 
efficiency of the tidal interaction is greatly enhanced. 

1.4. Inertial waves 

The study of inertial waves is said to have originated with Poincare (see Greenspan 
1968). Much attention has been devoted to the oscillations of an incompressible fiuid con- 
tained in a rotating spherical, spheroidal, or cylindrical container. When the Ekman number 
is small, the Coriolis force provides the only significant restoring force, away from viscous 
boundary layers, and the resulting inertial waves have remarkable mathematical and physi- 
cal properties. The low-frequency oscillations of an adiabatically stratified star or planet in 
uniform rotation are governed by a very similar problem (Papaloizou & Pringle 1981, 1982). 
In each case the inertial waves have frequencies (as observed in the rotating frame) no greater 
in magnitude than twice the angular velocity of the fluid, and have a rich spectrum that is 
dense or continuous in the absence of viscosity. 

In the simplest problems, such as that of a full sphere of incompressible fluid, an appar- 
ently complete set of inviscid inertial wave solutions can be obtained analytically (Greenspan 
1968; Zhang et al. 2001). However, the calculation of inertial waves in a spherical annu- 
lus, or in a star or planet with an arbitrary density stratification, requires two-dimensional 
numerical computations. Moreover, since the inviscid Poincare equation is spatially hyper- 
bolic in the relevant range of frequencies, the problem is ill-posed unless an explicit viscosity 
is included. Rieutord & Valdettaro (1997) calculated inertial waves in a spherical annu- 
lus, noting the behavior of the frequencies, damping rates, and eigenfunctions in the limit 
of small Ekman number. Their results were further elucidated by Rieutord, Georgeot, & 
Valdettaro (2001), who explained the important influence of the sohd core. In particular, 
the viscous eigenfunctions are concentrated on characteristic rays of the Poincare equation, 
which typically converge towards limit cycles as they reflect repeatedly from the inner and 
outer boundaries. A balance between the focusing of the wave energy along the converging 
rays and a lateral viscous diffusion sets the characteristic width of the eigenfunctions and 
determines the damping rates of the modes. 

Dintrans & Ouycd (2001) used similar methods to calculate inertial waves in a slowly 
rotating polytrope at more modest Ekman numbers, with potential application to Jovian 
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seismology. They considered a pure inertial wave problem, eliminating acoustic and gravity 
waves by adopting the anelastic approximation and using a neutrally buoyant model. Indeed, 
in contrast to either high-mass or solar-type stars in close binaries, the entire envelope of 
a gaseous giant planet is likely to be convective except for a shallow atmospheric layer. In 
contrast to the assumption adopted by loannou & Lindzen (1993b), efficient heat transport 
by convection prevents the propagation or excitation of g-mode oscillations. Inertial modes, 
however, can propagate throughout the planetary envelope. 

Jupiter's spin frequency is a significant fraction of its dynamical frequency and is more 
than twice the orbital frequencies of the Gafilcan satellites. Short-period extrasolar planets 
probably formed as rapid rotators similar to Jupiter and Saturn, although most of them 
may have already estabhshed a state of synchronous rotation. In each case the tidal forcing 
occurs in the frequency range of inertial waves, which will therefore constitute the natural 
and dominant response of the planet. Moreover, in this frequency range, the Coriolis cou- 
pling between the radial and angular motions is strong and the traditional approximation is 
inappropriate. Therefore two-dimensional computations are required. 



1.5. Plan of the paper 

In this paper, we revisit the issue of the excitation and dissipation of dynamical tides 
within gaseous giant planets. The response of the planet to low-frequency tidal forcing 
separates naturally into an equilibrium tide and a dynamical tide. In the convective regions 
of the planet the dynamical tide takes the form of inertial waves and we consider a turbulent 
viscosity associated with convective eddies to act on the tidal disturbance. If the planet 
has an outer radiative layer, outwardly propagating Hough waves are also excited and we 
consider them to be damped in the atmosphere. We therefore allow, in principle, for three 
avenues of tidal dissipation: viscous dissipation of the equilibrium tide, viscous dissipation 
of inertial waves, and emission of Hough waves. Our analysis permits the planet to rotate 
differentially, although we focus on the case of uniform rotation in our numerical calculations. 

In Section 2, we briefly recapitulate the different components of the tidal potential that 
are responsible for evolution towards synchronization of the planet's spin and circularization 
of its orbit. In Section 3, we carry out an analysis of low-frequency oscillation modes in 
a rotating planet. We formulate the perturbation equations governing both the convective 
and radiative regions of the planet's interior and analyze the matching conditions between 
these regions. In Section 4, we consider the tidally forced disturbances, including the equi- 
librium tide, the dynamical response in both regions, and the matching across the interface. 
We construct, in Section 5, two numerical schemes to obtain global solutions of the forced 
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perturbation equations. In Section 6 we present and discuss numerical results appropriate 
for the internal structure of gaseous giant planets. In Section 7, we compare our results 
with those obtained in previous investigations. We summarize our findings and discuss their 
implications in Section 8. 



2. COMPONENTS OF THE TIDAL POTENTIAL 

We consider two bodies in a mutual Keplerian orbit, and adopt a non-rotating coordinate 
system with origin at the center of body 1. In our study body 1 will be a giant planet and 
body 2 its host star. Let r be the position vector of an arbitrary point P in body 1, and let 
R{t) be the position vector of the center of mass of body 2. When evaluating the tidal force 
acting within body 1 it is usually adequate to treat body 2 as a point mass, i.e., to neglect 
its rotational and tidal deformation. The potential of body 2 experienced at the point P, 
including the indirect potential resulting from the acceleration of the origin of the coordinate 
system, is then 

GM2 GM2,^ , GM2 GM2 r„2 2 o/^ ^/ 



\R-r\ V / 2R 

The tidal potential is usually considered to be the non-trivial term of lowest order in r. 



* = ^[i?V-3(i^.rf]. (2) 



Without loss of generality, let body 2 orbit in the xy-plane, so that its Cartesian coor- 
dinates are 

i? = i?(cos(y9, sin V?, 0), (3) 
while the Cartesian coordinates of the point P are 

r — r (sin 9 cos 0, sin 9 sin 0, cos 9) , (4) 

where (r, 9, (f)) are spherical polar coordinates. Then 

* = ^^r" [l-3sin2^cos2(0-<^)] . (5) 



For a slightly eccentric orbit of semi-major axis a, eccentricity e, and mean motion n, 
and with a convenient choice for the longitude and time of pericenter passage, we have 



= a(l - ecosnt) -\- O(e^), ^p^nt-\- 2es\nnt -\- O(e^). 



(6) 
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The tidal potential may then be expanded in powers of e as 

* = ^(0) + ^(i) + O(e2), 

where 



Let 



GM2 2 



{2-3 sin^ ^ [1 + cos(2(/) - 2nt)\ } , 



(7) 
(8) 



4a3 

er^ {4 cos nt + sin^ [cos(20 — nt) — 7 cos(20 — 2>nt) — 6 cos nt] } . (9) 



Pr(cos^) 



(2£ + l)(£-m)!l'/' 



Pr(cos^), 



2(£ + m)! 

where < m < ^, denote an associated Legendre polynomial normalized such that 



sinOde = 1. 



(10) 



(11) 



The tidal potential correct to 0(e) can then be expanded in a series of rigidly rotating solid 
harmonics of second degree, 



^(0) + ^(1) ^ ll^Re Ajr'P^^ (cos e)e 



i{m,j<p—uijt) 



(12) 



where 



10 



mi = 0, uji — 0, Ai 

m2 = 2, a;2 = 2n, ^2 = - 

ms = 0, UJ3 — n, As — Be 

777.4 = 2, a;4 = n, A4 — e 

1715 — 2, UJ5 — 3n, 



1/2 



10 



1/2 



1/2 



20 



1/2 



1/2 



(13) 



Zahn (1977) gives an equivalent expansion carried to O(e^). 



When body 1 is a differentially rotating planet with angular velocity Q(r, 9), the effec- 
tive frequency of component j as experienced by the rotating fluid is the Doppler-shifted 
frequency cuj — cuj — rrijD,. Component 1 can be neglected because it is independent of time 
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and causes only a hydrostatic distortion of the planet but no tidal dissipation. The effective 
forcing frequencies of the remaining four components are 

a)2 = 2(n — il), cus — n, 0)4 = n — 2Q, 0)5 = 3n — 2Q. (14) 

This prescription assumes that the rotation axis of the planet is perpendicular to its or- 
bit. Otherwise an additional component of the tidal potential appears at first order in the 
obliquity, having ttiq — 1 and cuq — n, and therefore cuq — n — D, — 002/'^- 

In the special case of a uniformly rotating planet, the forcing frequency of component j 
lies in the spectrum of inertial waves if —2 < ujj/D, < 2. In this frequency range, the Coriolis 
coupling between the radial and angular motions is strong (Savonije et al. 1995). For any 
planet being spun down by tidal dissipation towards a synchronized state (i.e., Q > n > 0), 
all four effective forcing frequencies lie within the spectrum of inertial waves, as illustrated 
in Fig. 1. Once synchronism is achieved (i.e., Q = n), component 2 ceases to cause tidal 
dissipation but, if the orbit remains elliptical, the three eccentric components of the tidal 
potential continue to provide forcing at lj = ±f2, clearly within the spectrum of inertial 
waves. 

These considerations arc directly relevant to the short-period cxtrasolar planets. All the 
known satellites of Jupiter and Saturn, except the very distant, retrograde ones, also provide 
low-frequency tidal forcing within the spectrum of inertial waves. We therefore suggest that 
the Coriolis force has a dominant influence on tidal interactions in many cases of interest. 
Except in work by Savonije and co-workers on tides in uniformly rotating stars (Savonije et 
al. 1995; Savonije & Papaloizou 1997; Papaloizou & Savonije 1997; Witte & Savonije 1999a, 
1999b, 2001, 2002; Savonije & Witte 2002), the Coriohs force has been unjustifiably neglected 
owing to the mathematical, physical, and computational difficulties that it introduces. 



3. LOW-FREQUENCY OSCILLATION MODES OF A ROTATING 

PLANET 

3.1. Basic equations 

We consider a planet that is predominantly gaseous, but may contain a solid core. We 
initially consider the fluid regions to be ideal, satisfying Euler's equation of motion, 

^ = --Vp-V$, (15) 
JJt p 

the equation of mass conservation, 

^ = -pV.., (16) 
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and the adiabatic condition, 



1 Dlnp Dlnp 



0, 



7 Dt Dt 

where u is the velocity, p the density, p the pressure, 7 the adiabatic exponent, and 

D d „ 

the Lagrangian time-derivative. The gravitational potential satisfies Poisson's equation. 



(17) 



(18) 



(19) 



When expressed in spherical polar coordinates (r, 9, 0), these equations take the form 



with 



DUr 

'Dt 



r 



r 



Due ^ UrUe cos 9 



Dt 



Du^ 



Dt 



+ 



+ 



Dp 
Dt 



-P 



r 

1 d 



r r sin ^ 

UQUff, cos 9 
r sin 9 



Idp _ 9$ 
pdr 9r ' 

1 dp 19$ 
pr 89 r 89^ 
1 8p 18^ 



pr sin 9 8(f) r sin 9 8(1) ' 



8 



8r r sin 9 89 

IDlnp Dlnp 



U0 sin 9) + 



rsin^ 90 



7 Dt 



Dt 



0, 



89)'"r''s\n^9 8(l>'' 



— + u — 
8t ^ 8r 



ue 8 



+ 



U4, 8 



r 89 r sin 9 1 



(20) 

(21) 
(22) 
(23) 
(24) 
(25) 

(26) 



3.2. Basic state 

The unperturbed state is a steady, axisymmetric planet which may rotate differentially 
in its fluid regions. We introduce the small parameter e, being a characteristic value of the 
ratio of the spin frequency to the dynamical frequency (GMi/i?f)^/^ of the planet, which 
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has mass Mi and radius Ri . Centrifugal distortions of the planet are of fractional order 
and so we may pose the expansions 

$ = $Q(r) + e2$2(r,^) + 0(e^). (27) 
We neglect any meridional circulation and set Ur — ue — 0. The angular velocity is 

"""^ =0 = 60(1) (r,^). (28) 



rsin^ 



From equations (20) and (25) at order e*^, wc obtain the standard equations of structure 
for a spherical planet. We define the (inward) radial gravitational acceleration g[r) and the 
Brunt-Vaisala frequency N{r) in the usual way by 

= ^ ^ ]_dpo f d\n.pQ _ 1 (ilnpo \ 

dr ' Pq dr \ dr ^ dr J 

Equation (22) is satisfied identically. In practice the profile of differential rotation is deter- 
mined by the meridional circulation and by viscous, turbulent, or magnetic stresses that do 
not feature in the equations as written above. We will not require a solution of the basic 
state at order or higher. 

Our analysis is based on asymptotic expansions in the small parameter e^. Jupiter is 
a relatively rapid rotator, having Q'^Rf/GMi 0.08. This number increases to about 
0.14 in the case of Saturn, but is greatly reduced, to about 0.0016, if Jupiter is spun down 
to a period of three days. This last case may be taken as representative of short-period 
extrasolar planets once they achieve synchronous rotation. 



3.3. Linearized equations 

We consider the linearized response to a tidal potential component of the form 

Re [*(r,^)e*(™^-^*)] , (30) 

where m e Z is the azimuthal wavenumber and a; e R the forcing frequency. As in Section 2, 
we adopt the convention that m > 0. The induced Eulerian velocity perturbation is 

Re [u'{r, 9) e^^"^'^"'^*)] , (31) 
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and similarly for other perturbed quantities. The linearized equations then take the form 

- 2Q sin ^ < = + 1 f p'l? - p'l? V (32) 



where 



dr p2 \ dr dr 
^ , n , ^dW I f ,dp .dp\ 

/ 1 f ' ^ \ /r> 2 • 2 /i\ imW 



-icou'^ + — < + ^ (l^r^sin^e) = — , (34) 

^ rsmu \ or r ou I rsmu 



_ , .dp Un dp 



(rX) + ^^T^Ksm^) + 



r"^ dr ^ r sin 9^ r sin ^ 



(5r V dr J sin 9 dO \ dO ) sin ^ 



(35) 



\7P P/ V7 dr J r \7 9^ 5^ / ' 

r-"^ + IT— 7:^ sm^— - -— ^ = 47rGp', (37) 



(2; = a; — mVL (38) 
is the Doppler-shifted forcing frequency, and 

= ^ + + ^ (39) 
P 

is a modified pressure perturbation. In adiabatically stratified regions, where Vlnp = 

7Vlnp, the Eulerian perturbations of pressure and density are related by p' /p = ^p' / p, 
and the terms in parentheses in equations (32) and (33) cancel. 



3.4. Free modes of oscillation 

We first consider the free modes of oscillation of the planet, neglecting the tidal forcing. 
We treat the convective regions of the planet as being adiabatically stratified (Vlnp = 
7Vlnp), and the radiative regions as being stably stratified such that the Brunt-Vaisala 
frequency is comparable to the dynamical frequency. Where a radiative region adjoins a 
convective region, however, there is a transition region in which the Brunt-Vaisala frequency 
goes to zero. We obtain asymptotic solutions for low-frequency oscillations in each region 
separately and then match them together. For clarity, we use accented symbols, e.g., a, a, 
and a, to distinguish solution components in convective, radiative, and transition regions, 
respectively. 
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As explained in Section 2, we are interested in oscillation frequencies of order e, compa- 
rable to the angular velocity of the planet but small compared to the dynamical frequency. 
Accordingly, we set u) = ecc;(i) . To achieve such a low frequency in the presence of potentially 
large restoring forces, the perturbations are highly constrained. They should be almost in- 
compressible (or, more correctly, anelastic) in order to avoid a large acoustic restoring force. 
In radiative regions the perturbations should also involve predominantly horizontal motions 
in order to avoid a large buoyant restoring force. In fact, the relevant low- frequency modes 
are inertial waves in convective regions and what we call generalized Hough waves in radia- 
tive regions. In the solid core, if present, the rigidity is such that no low-frequency modes 
are permissible. 



3.5. Convective regions 

In convective regions we seek solutions of the unforced linearized equations of the form 

< ~ eu'^{r,e), 

p' ~ e'p'{r,e), 
p' ~ e'p'{r,0), 

~ e^$'(r,^). (40) 

The overall amplitude of the perturbations is arbitrary, but the scaling adopted here is 
relevant to the tidally forced problem considered in Section 4 below. The relative scalings of 
the perturbations are chosen to satisfy the linearized equations of Section 3.3, which reduce 
at leading order to 

-ia)(i)ti^ - 2fi(i) sin6'M^ = (41) 

1 dW 

-iuj (i)u'q- 20.(1) cos 9 ii'^ ^ (42) 

; 1 f , d u'o d \ ,^ 9 . 9 imW 

-icu^Du'^ + —— it'— + ^— {n(,)r^sm^e) = — , 43 

y ' V r sm \ or r oO J rsmO 

'r'poK) + — -7;^(«^sin^) + —4 = 0, (44) 



r^po dr ^ r sin 9 d9 r sin 9 

Po Po 



(45) 
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with 
and 



1 d f ^d^'\ 1 d f . M'\ m2$' ^ 

+ ^^^^^ ~ 2 • 2. = 46 

r"' ar \ ar / r"' sm oO \ oO J sm^ 

'^(i) = ^(1) - "^Q(i), (47) 

= ^ + (48) 

Po 

Equations (41)-(44) are decoupled from those that follow, and, when combined with appro- 
priate boundary conditions, define the inertial wave problem. We note that equation (44) is 
equivalent to the constraint V • (pow') = of the anelastic approximation (Ogura & PhiUips 
1962), which eliminates acoustic waves from the problem. Self-gravitation also plays no role 
for inertial waves, because simply combines with p' in the combination W. 



3.6. Radiative regions 

In radiative regions of the planet, the oscillation frequency is small compared to the 
Brunt-Vaisala frequency, implying that the motions should be predominantly horizontal in 
addition to being anelastic. These constraints are satisfied when the radial wavelength is 
short compared to the radius of the planet, and the perturbations then have a WKB form 
in the radial direction. Accordingly, we seek solutions of the unforced linearized equations 
of the form 

< ~ e''u',{r,e)E{r), 
u', ~ eu'e{r,e)E{r), 
u'^ ~ eu'^{r,9)E{r), 
p' ~ ep'{r,e)E{r), 
p' ~ e^p'{r,e)E{r), 

~ e^$'(r,^) + e^l>'(r,^)^(r), (49) 

with 

(50) 

The exponent s < remains arbitrary until the solutions in adjacent convective and radiative 
regions can be connected. We show in Section 3.7 below that this procedure yields s — —1/6. 
Note that the term e^$'(r, 9), which does not have a WKB form, is the smooth continuation 
of the gravitational potential perturbation associated with the inertial waves in convective 
regions. 



E{r) 



e^'exp 



k{r) dr 
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The linearized equations of Section 3.3 then reduce at leading order to 

9P' 



= -ikW - 



Po 



— ia)(i)it^ — 2Q(i) cos 9u'^ — — 



IdW 



sin 9 do 
1 d 

iku'^ + — - — {ug sin 9) + 
r sm 6 oO 



r 89 

imW 
r sin ^ ' 

= 0, 



r sin 9 



with 
and 



. P' 
icu/i) \ ii„ = 0, 

Po 9 



Po 



(51) 

(52) 

(53) 

(54) 

(55) 
(56) 

(57) 
(58) 



Since $' appears only in Poisson's equation, self-gravitation is irrelevant for waves of such 
short radial wavelength. The remaining equations may be combined into the form 



U2 2 



where £ is a linear operator defined by 



d 2mQ(^i)Cos9 
sin 9 \d9 a)(i) sin^ 



7V2 



sin 6* / d 2mQ(i) cos 9 



D \d9 



tD(i) sin 9 



w 



+ 



a)^-^^ sin^ 9 ' 



with 



2 .o2 _2.. OO '^^Wcos^sin^. 



D = cu(^) - 4f]^i) cos" 9 - 20(1) 



Let {Xi} and {wi{9)} be the eigenvalues and eigenfunctions of C, satisfying 



(59) 



(60) 



(61) 



(62) 



If the angular velocity ^2(1) is independent of r, the eigensolutions will be also. For a general 
pattern of differential rotation, however, the eigensolutions depend parametrically on r. In 



- 21 - 



either case, each eigensolution with A > corresponds to a freely propagating wave solution 
(which can be considered a rotationally modified g mode) satisfying the dispersion relation 



For a given wave frequency u and azimuthal wavenumber m, this relation determines the 
radial wavenumber k at each r. Evanescent eigensolutions with A < are also possible in a 
rotating planet. 

As given above, the solution represents a traveling wave with phase velocity e^uj{i)/k. A 
standing wave is formed from a superposition of two such solutions with opposite choices of 
the sign of k. Depending on how the wave is excited, and on whether it is allowed to refiect 
from the boundaries of the radiative region, the relevant solution may be a traveling wave, a 
standing wave or, in the case of partial reflection, an intermediate solution. In the case of a 
traveling wave the sign of k should be chosen with due regard to the direction of the radial 
group velocity. 

In the absence of dissipation, the waves carry conserved fluxes of energy and angular 
momentum, owing to the invariance of the basic state with respect to t and (j). These 
conservation relations can be derived from the theory of 'wave action' using a Lagrangian 
formulation of ideal fluid dynamics (e.g., Hayes 1970). For a traveling wave, the average flux 
of energy wave action through a spherical surface can be shown to be 



and the flux of angular momentum wave action differs only by a factor m/uj(i). Outwardly 
propagating waves, having positive radial group velocity, are those with sgn(a)(i)/A;) < 0. 
This property occurs because the frequency of g modes decreases with increasing radial 
wavenumber. In the absence of dissipation the flux F should be independent of radius. 

We note some properties of the operator C Singular points of the equation Cw — Xw 
occur at the poles and wherever a)(i) = ('corotation resonance') or D = ('Lindblad 
resonance'), although the latter is only an apparent singularity. Close to the north pole, two 
independent solutions are w ~ ^"^ and w ~ 9~'^ (or w ~ 1 and w ~ In^ in the case m — O). 
The condition that w be bounded selects the regular solution. This, and a similar condition 
at the south pole, provide the two boundary conditions for the eigenvalue problem. 



F 




(64) 



(65) 



In a differentially rotating planet, a corotation resonance can occur, in principle, on any 
surface of constant angular velocity. As noted by Goldreich & Nicholson (1989), this feature 
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is expected to present an absorbing barrier to propagating waves. In this paper we disregard 
the possibihty of corotation resonances, but their influence merits further investigation, which 
we defer to future work. 



The operator L is self-adjoint in the sense that 

/»7r r /»7r 

Jo Uo 



de 



(66) 



for any two functions u{9) and v{9) satisfying the regularity conditions at the poles. The 
eigenvalues A are therefore real and the eigenfunctions corresponding to distinct eigenvalues 
are orthogonal in the sense that 



w*Wj sin 9d9 — 0, i ^ j. 



Furthermore, if w{9) is a function satisfying the regularity conditions. 



(67) 



/ w* (Cw) sin 9 d9 ^ 
Jo JO 



1 

D 



d 2mfi(i)Cos^^ 
89 a)(i) sin 9 



w 



+ 



21 12 

m \w\ 



u!7, N sin 9 



sm9d9. 



(68) 



This equation shows that the eigenvalues of C are all positive when D > throughout 
< ^ < TT. If D changes sign in < ^ < tt, there will be sequences of positive and negative 
eigenvalues. 

Later, we will need to know under what conditions the operator C can be inverted. 
Provided that the eigenfunctions form a complete set within the relevant function space, the 
inversion is possible unless C has a null eigenvalue A = 0, in which case the inversion can be 
performed only within the subspace of functions orthogonal to the null eigenfunction. In the 
case 171 = 0, the only possible null eigenfunction is Wq = constant, which exists for all values 
of u;(i). In the case m ^ 0, this solution does not occur, but there may nevertheless be null 
eigensolutions for isolated values of uj(^iy 

In the case of a non-rotating planet, £ reduces to a multiple of the horizontal Laplacian 
operator and the eigensolutions are 



A,; = 



n{n + l) 



UJ 



W^{9) 



pr(cos« 



n 



m + i, 



i > 0. 



(69) 



(1) 



In the case of a uniformly rotating planet, or one in which the angular velocity is independent 
oi 9, jC reduces to Laplace's tidal operator and the eigenfunctions are the Hough functions 
(e.g.. Chapman & Lindzen 1970). We note that, in the case m = 0, Laplace's tidal operator 
has null eigenfunctions when = — 2mQ(i)/n(n + 1) for integers n > m. These solutions 



-23- 



are known variously as planetary waves, Rossby waves, r modes, or toroidal modes, and 
correspond to global modes (hence A; — > 0) of a uniformly rotating radiative region (e.g., 
Papaloizou & Pringle 1978). 

For a general pattern of differential rotation in which Q depends on 9, we call the 
eigensolutions 'generalized Hough modes'. The general solution in radiative regions involves 
a linear combination of generalized Hough waves. 



3.7. Matching between convective and radiative regions 

In planets and stars, convective and radiative regions can be disposed in a variety of 
configurations, and planets may also contain solid cores. We consider in detail the case 
in which a single radiative envelope exists exterior to a single convective region, possibly 
containing a solid core. This is relevant to short-period extrasolar planets and also to high- 
mass stars, although in the former case the radiative envelope is shallow. The same principles 
can be used to construct solutions for alternative arrangements of convective and radiative 
regions. 

Near the boundary r — between convective and radiative regions, we assume that N"^ 
is zero on the convective side and rises linearly from zero on the radiative side, so that 



V 



dr 



> 0. (70) 

r=rb+ 



As r — > Tb from above, the wavelengths of all the generalized Hough modes tend to infinity 
because k oc N, and the approximations adopted in Section 3.6 break down. The situation is 
related to the problem of Lindblad resonances in differentially rotating discs (e.g., Goldreich 
& Tremaine 1978; Lin & Papaloizou 1979a; Lubow & Ogilvie 1998). By analogy with that 
problem, the characteristic radial extent of the transition region can be identified as being 
of order e^/^. We therefore write 

r = Tb + e^/^x, (71) 
where x is an inner variable, of order unity within the transition region. We then have 

= e^/^Vx + 0{e''/^). (72) 

We pose the following expansions for the perturbations within the transition region: 

< ~ e^/%{x,e), 
u'g ~ eug(x,e), 
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^0 ~ eu'^(x,e), 

p' ~ e'/'p'ix,e), 

~ e2|''(rb,^) + e«/^l''(x,^). (73) 

Note again that the term e^<l>'(rb,^^) is the smooth gravitational potential perturbation as- 
sociated with the inertial waves in convective regions, evaluated at the boundary. 

The linearized equations of Section 3.3 then reduce to 

^ - ^. (74) 

OX po 

1 dW 

-ia)(i)Mg - 2f](i) cos6'm^ = (75) 
/ u'f, d ,^ . 9 imW 

1 |«.in.) + ^ = 0. (77) 

OX r sm d o6 r sm 

o' Vx 

icOm— + < = 0, (78) 

Po 

= 47rG'p', (79) 



with 



dx^ 

W = ^ + ^'. (80) 

Po 

In these equations all coefficients are to be considered independent of x, except where it 
appears explicitly, and are to be evaluated at r = rb. Elimination of perturbations in favor 
of W now yields 

r2 d fldW 



~^ V dx \ X dx J ^' ^ 
where, again, all coefficients are to be evaluated at r = rb. This is the equivalent of equation 
(59) in the transition region. 

Equation (81) can be solved by separation of variables. Let Hi be the real solution of 

4 = (82) 
Then the general solution of equation (81), regular at the poles, is 

W = Y^ = J2^i [^i'iKix) + Sii Bi'{Kix)] Wi{e), (83) 
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where {oj} and {sj} are arbitrary coefficients, and Ai' and Bi' denote tfie derivatives of tfie 
Airy functions of tlie first and second kinds. Tliis solution must be matcfied to tfiat in tfie 
radiative region at a; ^ 1, and to tlie solution in the convective region at a; = 0. 

For generalized Hough modes with A > 0, k is negative and the wave propagates in 
the radiative region x > 0. If the boundary condition at the outer surface of the radiative 
region is such that the waves are not reflected, the relevant solution is a traveling wave with 
outwardly directed group velocity. This is achieved by setting either Sj = ±1, because the 
asymptotic form of the Airy functions is such that, for < 0, 



Ai'{Kx) ±iBi'{Kx) 



-TT ^/^(-ra)^/^ exp <^ 



li-^xf' + l 



(84) 



as X — > oo. This matches on to the WKB form of the solution in the radiative region, because 
the wavenumber can be identified as 



-2/3 



dx 



-i-Kxf/^ + ? 



(85) 



in agreement with the WKB wavenumber e given by equation (63). In fact the relevant 
solution has Si = sgn(a)(i)). 

For modes with A < 0, k is positive and the wave is evanescent in the radiative region. 
The relevant solution is then Si = 0, because Ai'(fi;a;) decays exponentially as kx — > oo, while 
Bi'(/ta:) grows exponentially. 

To match the solution in the transition region to that in the convective region at x = 0, 
we first compare the ordering schemes (40) and (73), and find that u'^. and p' arc formally 
mismatched in order. Now equation (78) implies that p' ^ as x ^ 0, which resolves the 
mismatch in p'. The matching of u'j. implies that the effective boundary condition for inertial 
waves in the convective region is that 



;(rb,^) = 0. 



(86) 



Matching of the other variables is achieved by requiring that W be continuous at r = rb, 
i.e., 

W(r^,e)^W(0,e). (87) 

Note that, if W is continuous across the spherical boundary, u'q and u'^ are also continuous 
because they are related to one another only by horizontal derivatives, and in exactly the 
same way on each side of the boundary. The matching condition 

^(^b, 0)=J2(^^ [Ai'(O) + s,t Bi'(O)] w0) 
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determines the amplitude coefficients according to tfie ortfiogonal projection 

Jq w*W{rk, 6) sine de 
" [Ai'{0)+s^iBi'{0)]£\wi\^ smOdO' 



(89) 



Tfie pliysical interpretation of tfie matcliing conditions is tfiat tfie transition region 
appears as a rigid boundary for inertial waves in tfie convective region because ttie steeply 
rising entropy gradient suppresses radial motions. The generalized Hough waves in the 
transition and radiative regions are excited by the modified pressure perturbation W acting 
at the boundary r — ry,. 

The outward energy flux in generalized Hough waves is found by evaluating equation 
(65) in the transition region at large x, with the result 

2\ 1/3 



Poc.a)sgn(c.a))^^J Y.\ r\^^\2,^ed9 



-(90) 



where the sum includes only those modes for which A > 0, and all quantities are to be 
evaluated at r = rb- Comparison with equation (65) indicates that s = — 1/6, as anticipated 
above. A more detailed matching between the transition region and the radiative region is not 
necessary for our purposes, but the constancy of the energy flux could be used to determine 
how the amplitudes of the Hough modes vary with radius throughout the radiative region. 

It is likely that the waves will become nonlinear and be damped in the outer layers of 
the planet. In this case the emission of generalized Hough waves at the boundary between 
convective and radiative regions provides a dissipation mechanism for the inertial waves. 
The rate of energy loss in Hough modes is formally of order e^^/'^, while the energy content 
of the inertial modes is of order e^. This route of dissipation must be compared with viscous 
or turbulent damping of the inertial waves in the convective region itself. 



3.8. Viscosity and inertial waves 

The procedure suggested by the foregoing analysis is to solve the inertial wave problem 
defined by equations (41)-(44) in the convective region, using the rigid boundary condition 
(86) at the outer surface. If the planet has a solid core, a rigid boundary condition is also 
appropriate there. In the absence of a core, regularity conditions apply at r = 0. The 



amplitudes of the generalized Hough modes in the radiative region can then be computed 
from equation (89). 

In a uniformly rotating planet, equations (41)-(44) admit a dense or continuous spec- 
trum of inertial wave solutions (e.g., Papaloizou & Pringle 1982) occupying the frequency 
interval — 2|fi(i)| < a)(i) < 2|Q(i)|. In a differentially rotating planet, a similar phenomenon 
occurs, although the extent of the spectrum is slightly different (e.g., Lin, Papaloizou, & Kley 
1993). This property is associated with the fact that the system of equations is spatially 
hyperbolic, rather than elliptic, for wave frequencies within the dense or continuous spec- 
trum. For such frequencies, the mathematical problem requires the imposition of Cauchy 
boundary conditions on an open surface, while physical considerations lead to rigid (non- 
Cauchy) boundary conditions on a closed surface. The inertial wave problem is therefore 
mathematically ill-posed in the absence of viscosity. 

When studying the tidally forced problem, we are interested in calculating the rate of 
dissipation of energy, or, equivalently, the tidal torque exerted on the planet. For effective 
forcing frequencies within the range of the dense or continuous spectrum of inertial waves, the 
possibility exists of resonantly exciting an infinite number of modes. We therefore introduce 
a small kinematic viscosity u in convective regions, partly to regularize the problem and 
partly in an effort to model the effects of turbulent convection on the wave motion. In 
the presence of viscosity, the inertial wave problem is mathematically well posed and the 
spectrum is discrete. As the viscosity tends to zero the spectrum becomes increasingly dense 
(e.g., Dintrans & Ouyed 2001) and the probability of resonant excitation of wave modes 
increases. It is therefore possible to envisage a situation in which the dissipation rate would 
not vanish linearly with i/ as — > 0. Indeed, it might even be independent of as — > 0, as 
occurs in hydrodynamic turbulence. 

In the presence of viscosity it is necessary to introduce additional boundary conditions 
for the inertial wave problem. If the viscosity is restricted to the convective region r < rb, 
then continuity of the velocity and of the stress requires that the rO- and r0-components of 
the viscous stress tensor vanish at r = rb. The boundary conditions at r = rb are therefore 

or \ r J or \ r J 
If the planet has a solid core of radius rc, the no-slip boundary conditions 

u', = u', = = (92) 

are appropriate at r = rc. 
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4. TIDALLY FORCED OSCILLATIONS 



4.1. The equilibrium tide 



The response of a body to tidal forcing is traditionally calculated using the theory 
of the equilibrium tide, in which an initially spherical body adjusts hydrostatically to the 
varying tidal potential (e.g., Darwin 1880). This approach is usually justified on the basis 
that the tidal forcing frequency is much less than the dynamical frequency of the body. By 
introducing an arbitrary phase lag into the response, usually parametrized by the quality 
factor Q, the theory allows for a tidal torque to be exerted and energy to be dissipated. The 
efficiency of known tidal interactions in the solar system places constraints on the Q-values 
of the bodies involved (Goldreich & Soter 1966). 

Such a theory may offer a reasonable description of the tidal response of solid bodies, in 
which the dissipation factor 1/Q reflects the imperfect elasticity of the medium. However, 
it is quite inappropriate for fluid bodies, in which the dissipation is of a different nature 
and, moreover, a wavelike response, the dynamical tide, occurs at low forcing frequencies 
(Cowhng 1941; Zahn 1966). Although, in principle, tidal dissipation in a fluid body can 
still be parametrized by Q, the required Q- value will depend on the tidal forcing frequency 
and on the properties of the body in such a complicated way that the parametrization is of 
questionable value. 

In modern studies, the response of a fluid body to tidal forcing is conveniently separated 
into the equilibrium tide, which represents a large-scale, quasi-hydrostatic distortion of the 
body (without introducing a phase lag), and the dynamical tide, which usually constitutes 

a relatively small, mostly wavelike correction that is critically important because of its con- 
tribution to the dissipation rate. The equilibrium tide associated with a tidal potential ^ is 
a Lagrangian displacement field ^ having the properties (e.g., Goldreich & Nicholson 1989) 



(^(e) + 



v-^ = o. 
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9 



The Eulerian perturbations associated with the equilibrium tide are 



■(e) = e [-^'^(i)^ - rsin6'(^ • V^^i))e^] , 



(94) 





(95) 



where ^'^^^ and p'^^^ are related by 



V%)=47rGp;,). 



(96) 
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Note that V • ^ = implies V • u'^g^ = 0, each implying that the equilibrium tide is incom- 
pressible.^ 

We consider a single component of the tidal potential in the form of a spherical harmonic 
of degree i > m, 

^if ^■^{r)Pp {cos 9). (97) 

As noted in Section 2, the more important components of the tidal potential have £ — 2 and 
m = 0, 1, or 2. The solution for the equilibrium tide is of the form 

sm u 

= |.;(r)Pr(cos^), (98) 
where the radial functions satisfy the ordinary differential equations 

9 

^^(^^«-^4 = o, COO) 

P'. = (101) 
A = -(/f. (102) 



dr \ dr I r 



The boundary conditions are that the solution should be regular at the center and at the 
surface of the planet. At the surface, l»g should also match on to a decaying sohd harmonic 
oc r''^^'^^\ with continuity in the function and its first radial derivative. 

Note that the equilibrium tide is not affected by rotation, and involves only a single 
spherical harmonic. This is not true of the dynamical tide. 



^In expressions such as V • ^ and V^$|g^ the azimuthal part of the operator is intended to be included, 
even though a Fourier analysis in ^ has been carried out. 
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In both convective and radiative regions, the perturbations p'^^^, p'^^^, and associated 
with the equihbrium tide constitute the dominant part of the response to tidal forcing, and 
it is convenient to subtract them out. In radiative regions it is also convenient to subtract 
out the radial velocity perturbation u'^^^y which dominates over the wavelike part. 



4.2. Convective regions 
In convective regions the forced solution is of the form 















rsj 








rsj 




p' 


- P[e) 


rsj 




p' 


-P'ie) 


rsj 








rv-/ 






w 







(104) 



The hnearized equations of Section 3.3 reduce to 



dW 

-iu(^i)u'^ - 2Q(i) sin 6u'^^ — —, (105) 

-iLo^i)Ug - 2fi(i) cos 61 -u^ = ^"^^^^ 

-lu^Du'^ + — u'— + ^— ^(Dr'sm^e = — , 107 

(^VoO + — -^T^Ksm^) + — = -^a;(l)e.-T^^. (108) 



r^Po dr ^ r sin 9 dO r sin 6* dr 

Equations (105)-(108) are identical to (41)-(44), except for what is effectively a forcing 
term on the right-hand side of equation (108). The unusual positioning of the forcing term 
in the equation of mass conservation reflects the fact that the dynamical tide is not forced 
directly by the tidal potential. Indeed, the tidal forcing is balanced in the equation of motion 
by the pressure and potential perturbations associated with the equilibrium tide. However, 
the time-dependence of the associated density perturbation provides an unbalanced term in 
the equation of mass conservation, which appears as the driving term for u' in the equations 
as written above. In fact, the equilibrium tide tt'^g^ satisfies equation (108) exactly, but fails 
to satisfy the equation of motion because of its contribution to the inertial terms on the 
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left-hand side. Therefore the dynamical tide is forced very indirectly, through these inertial 
terms. 

Note that the dynamical tide u' — n'^g^ in convective regions is anelastic, hke the free 
oscillations studied in Section 3.5, while the equilibrium tide is incompressible. The total tide 
is therefore neither anelastic nor incompressible, but satisfies V • {pou') — w'^^^ • Vpo, which 
is equivalent to equation (108). In a similar way, although self-gravitation is unimportant 
for the dynamical tide, which therefore satisfies the Cowling approximation, it cannot be 
neglected for the equilibrium tide. 



4.3. Radiative regions 

In radiative regions the forced solution is of the form 



<- 






e%{r,e)E{r), 
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- P'ie) 




ep\r,e)E{r), 
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-P'ie) 




e^p'{r,e)E{r), 








e^^'{r,e) +e^^'{r,e)E{r), 




w 




e^W{r, e)E{r), 



with E{r) as in Section 3.6. The linearized equations then reduce exactly to equations 
(51)-(58). This implies that the dynamical tide is not directly forced in radiative regions. 



4.4. Matching between convective and radiative regions 

In a transition region of the kind considered in Section 3.7 the forced solution is of the 

form 













~ eu'g{x,9), 








p' 


- P[e) 
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-P[e) 


~ e'p'{x,0), 




- 
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(110) 



The linearized equations reduce to 



dW 
dx 



9P_ 
Po 



— ia)(i)M^ — 2f](i) cos9u'^ — — 



IdW 



r de ' 



r dr 

19... , _ dul 1 



+ 



r sin 9 89 

lijJiX) \ = 

Po g 



sm9d9 

d ,_, . imu', 
Uq sm 9) + 



rsin^ 



imW 
r sin ^ ' 

= 0, 



0, 



dx"^ 



with 



Po 



47rGp', 



P(e)P2 



p§ 



(Ill) 

(112) 
(113) 
(114) 
(115) 

(116) 
(117) 



As in Section 3.7, all coefficients in equations (111)-(115) are to be considered independent 
of X, except where it appears explicitly, and are to be evaluated at r = rb- The equations 
now contain forcing terms (those involving ^r) that are also independent of To deal with 

these, it is convenient to split the sohition into a 'particular solution', which is independent 
of X and satisfies equations (111)-(115) as written, and a 'complementary function', which 
depends on x and satisfies equations (111)-(115) with the forcing terms omitted. Thus 



W{x, 9) = W^\9) + W^''\x, 9), 



(118) 



and similarly for u'g and u'^, while u'^ and p' do not appear in the particular solution. The 
ordinary differential equations defining the particular solution can be combined, if desired, 
in the form 



where 



1 d 



d 



2a;(i)Q(i) sin^ ^cos6' d 



Cj{\) dr a)(i) sin 9 89 

is a forcing combination independent of x. 



D 



(119) 



-^e.^(r 1^(1)) 



(120) 
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For the complementary function, elimination of perturbations in favor of W^'^^ again 
yields 

V dx \x dx J (-'-^-'-) 

with the solution 

W^W ^Y^ai [M{Kix) + sn Bi'imx)] Wi{e), (122) 

i 

where Sj = sgn(c(;(i)) for propagating modes with Aj > 0, and Sj = for evanescent modes 
with Aj < 0. The wavelike part of the solution dominates the particular solution at large x. 

Matching W to the solution in the convective region at a; = determines the amplitude 
coefficients according to 

!o < [w-(^b, 0) ~ w^^)] sin e dd ^^^^^ 



' [Ai'(O) + Sii Bi'(O)] |wi |2 sin 9 d9 
The outward energy flux in generalized Hough waves is then 



Poa;a)Sgn(a;(,))^-J ^ A, T^^p" 



sin 6 d9 

(124) 

where, again, the sum includes only those modes for which A > 0, and all quantities are to 
be evaluated at r = rb- 

The physical interpretation of this analysis is that the particular solution W^'^'' repre- 
sents the bulk response of the transition region to residual tidal forcing (i.e., the failure of 
the equilibrium tide to provide an exact solution to the problem), while the wavelike comple- 
mentary function is its response to the pressure of the inertial waves acting at the boundary 
r = Tb, as occurs in the case of free oscillations (Section 3.7). 

As discussed in Section 3.6, the inversion of the operator C to solve equation (119) for 
\^ip) requires careful interpretation in cases in which C has a null eigenvalue. In the case 
m = 0, it is easy to show that is orthogonal to the null eigenfunction wq — constant. The 
solution W^^ therefore exists but is non-unique, as any multiple of wq can be added to it. 
This ambiguity affects only the amplitude of this non-propagating mode, however, and not 
the energy flux in propagating Hough modes. In the case m 7^ 0, if £ has a null eigenvalue 
for an isolated value of oj^i), and if has a non-zero projection on to the null eigenfunction, 
the corresponding mode will be resonantly excited and the present analysis breaks down 
(formally, W^'^'^ 00). This phenomenon, called 'toroidal mode resonance' by Papaloizou 
& Savonije (1997), is discussed further in Section 5. 
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The matching of u'^ implies that the boundary condition for tidally forced inertial waves 
in the convective region is 

<-e"X(e) = at r = rb, (125) 

The transition region therefore appears as a rigid boundary for the dynamical tide in the 
convective region, while the equilibrium tide continues smoothly across the boundary. In 
the presence of an effective viscosity confined to the convective region, the equivalents of the 
boundary conditions (91) and (92) are then 

<-^-Xw = |(^)=|(^)=0 at r = r„ (126) 

K - e"^<(e) ^u'q- e~^4(g) = u'^ - e~\(^^) = at r = Tc- (127) 

The assumption underlying the boundary condition at r = Tc is that the core is distorted 
exactly according to the equilibrium tide, and so the dynamical tide satisfies a no-slip con- 
dition. This approximation neglects the influence of the rigidity of the core in modifying its 
tidal distortion, which is however only a small error for typical cores of several Earth masses. 



5. NUMERICAL ANALYSIS FOR UNIFORM OR 'SHELLULAR' 

ROTATION 

5.1. Planet model and equilibrium tide 

In the remainder of this paper we construct detailed numerical solutions for a simple 
model of a planet. Where possible, we omit the ordering parameter e and the ordering 
subscripts that were needed in deriving the equations in previous sections. 

Our numerical method is suited to the case of 'shellular' rotation in which the angular 
velocity fl{r) is independent of latitude (Zahn 1992). Of course, this prescription includes 
the simple case of uniform rotation, which will be the main focus of our initial investigation. 
We consider the planet to be fully convective except for a small solid core, r < r^, and a 
shallow radiative envelope, rb < r < We model the planet as a polytrope of index 1, so 
that the density is 

(71 Ml \ sin kr 
w)— ■ 

where k is related to the radius Ri of the planet by kRi — n. The gravitational acceleration 
is 

g — ( sin kr — kr cos kr). (129) 
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While it is common practice to model fully convective giant planets as polytropes of index 1, 
we neglect any modifications of this structure within the convective region associated with 
the presence of a core and a radiative envelope. 

The tidal potential is of the form 

* = -^Ar^, 130 

where A is a dimensionless constant, as in Section 2. The equations for the equilibrium tide 
(Section 4.1) imply that 

v2($; + ^) = ^^(*; + (131) 

For a polytrope of index 1 this simplifies to the Helmholtz equation 

(V2 + fc2)($; + ^) = 0, (132) 
with the solution (regular at r = 0) 

^'^ + ^^^BRrr-"'h/2{kTl (133) 

where i? is a dimensionless constant to be determined. At r = $g must match smoothly 
on to a decaying solid harmonic oc r~'^^^^\ and therefore 



d\nr 

This condition determines the ratio B/A — 5/-\/2. 
The detailed solution for the equilibrium tide is 



= -3 at r = Ri. (134) 



= |M [(3 - A;V) sinA;r - 3A;rcos A;r] - A , (135) 

3 r A (3 — k'^r^) sin kr — 3kr cos kr 



Ml J \ a J kP'T 



sin kr — kr cos kr 



(136) 



We note that, if a less idealized density profile is adopted, the equilibrium tide must be 
determined by solving the ordinary differential equation (131) numerically. 
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5.2. Projection of the equations on to spherical harmonics 



The basic equations of Section 4.2 governing tidally forced inertial waves in the convec- 
tive region are 



. dW Fr 

—luJUr — zilsmOus = — 1 , 

or p 

^ ^ IdW Fg 

—iujue — 2il cos Ous = — — H , 

r d9 p 



—iCjUfj) + yiVt + — — I sin 9ur + 2Vt cos 9 uq 
I d , . 1 d 



dlnr . 

[r pUr) H r ' 



imW F^ 
rsm9 p ' 



{u0sm9) + 



dlnp 



(138) 
(139) 
(140) 
(141) 



r'^pdr^ ' ' ' rsva.9 89^ " ' ' rsva.9 ^' dr ^ 

where we spcciahze to the case of shellular rotation and omit all unnecessary subscripts and 
superscripts. Included here, as explained in Section 3.8, is a viscous force with components 



Fg 
F^ 



sin 9 
1 

sin 9 
1 

sin^ 9 



d d 

— (r^ sin 9 T„) + (?^ sin 9 T^e) + imrT^, 
or o9 

d d 

— (r^ sin ^^re) + {r^ sin 9 Tgg) + imr'^Tg, 



r r 
Tj,j, cot 9 



d 

— {r^ sin' 9 Tr^) + — (r' sin" 9Tg^)+ imr^ sin 9 



where 



T = 2//S + //b(V-u)l 

is the Navier-Stokes viscous stress tensor, with 

S = ^[Vw+(Vu)T -^(V-u)l 



(142) 



(143) 



(144) 



being the traceless shear tensor, p{r) the (dynamic) shear viscosity and p\,{r) the bulk 
viscosity. The components of T are 



T = 



Trd — 



T, 



rcf) 



Tge 



2p 
1^ 



dUr 



+ (/Xb - \p)V ■ u, 



dr 

d fU0\ 1 duj. 
dr \ r 

4 (-) + 

or \ r / 



r 89 
imur 
r sin 9 



^ , 1 8ug Ur \ , 9 N „ 

2/^(-9/ + v) + (/^b-iA^)V-w, 
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Te^ — n 



sin^ d 



r 86 \sm.6 

imUff, Ur 



with 



rsinO 
1 d 



rsinO 
uq cot 6 



d 



r sm ^ r sm 



(145) 



(146) 



Following the basic idea of Zahn (1966), we decompose the velocity into 'spheroidal' 
and 'toroidal' parts and then project all variables on to spherical harmonics such that 



ue 
W 



J]a„(r)Pr(cos 



im d' 



= 5^2rf„(r)P™(cos 



Pr(cos^), 
P-(cos^), 



(147) 



where the summations are over integers n > m. Even if the tidal forcing and equilibrium 
tide involve only one spherical harmonic, the Coriolis force couples spherical harmonics of 
different degrees. Note that c„ is the amphtude of the 'toroidal' part of the velocity, which 
is not present in the equilibrium tide but is induced by rotation and is incompressible. 

To project equations (138)-(141) on to spherical harmonics, we require the recurrence 
relations 

cos^P-(cos^) = g„+iP- i(cos^) + g„P- i(cos^), (148) 



sin^^4c'(cos( 
d9 



-- ng„+iP„™ i(cos^) - (n + l)g„P;;ii(cos ^^), 
for normalized associated Legendre polynomials, where we define the coefficients 



(In 



n 



m 



All? - 1 



1/2 



(149) 



(150) 



and it is to be understood that P^{cos9) — ior n < m. To deal with the angular 
components of the equation of motion, equations (139) and (140), we take combinations 
that correspond to 'divergence' and 'curl': 

1 d 



r sin 9 89 

im 



[sin 9{—i(jjue — 2Q cos 9 

dVt 



rsn\9 
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—iCju^ + ( 20 + 
1 d 



sin 9 89 
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dlnr ^ 

8W\ m^W 



89 



sin 9ur + 2fl cos 9 uq 
1 8 [Ff 



sm^9 



+ 



+ 



im 



rsm9 89 \ p J rsm9 \ p 



, (151) 
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1 d 
r sin 9 39 

im 



sin 9 



-iujua, + 2iZ + — — 
\ a m r 



sin 9ur + 20: cos 9 ug 



rsm9 



[—icuug — 2Q cos 9 u^) 



1 d fF, 



im 



rsm9 d9 \ p J rsm9\p 



(152) 



Equation (138) then gives 



-iLJUn - 2imVtrhn + 2Vlr [{n - l)g„c„_i - (n + 2)g„+iCn+i] 

^ da 
2/x 



dr p dr 



-n{n+l) 



P 



p dr \ r 
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dr V 



(153) 



(154) 



Equations (151) and (152) give 
—n{n + l)iLjbr^ 



im / ^ dfl 
I 20 + 



r V 



dlnr 



a„ — 2imQ,br, 



idr, 
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4 """rt 

dr 



-2[n(n+ 1) - 1] 



pK 



(155) 



and 



-n(n + l)iLOCn + - (20 + -^^] [{n + l)g„a„_i - nqn+ian+i] 
r \ d m r / 

-2Q [(n -l){n+ l)qnbn-i + n{n + 2)g„+i6„+i] - 2imQ.Cn 

n(n + l) Jl d / 4dc„\ oi/^c„\ 



(156) 



respectively. Finally, equation (141) gives 

1 d 



r'^p dr 

where 5 is the Kronecker delta function. 



r pan) — n{n + l)bn = —iu^, 



dlnp 
dr 



(157) 



The couplings are such that, if ^ = m = 2, the quantities a„, bn, and d^ are non-zero 
only for even values n > 2, while c„ is non-zero only for odd values n > 3. li i — 2 and 
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m — 1, the quantity ci is also non-zero. Finally, ii £ — 2 and m — 0, the quantity do is 
also non-zero, but Qq and bo can be taken to vanish, and equations (155) and (157) are not 
required in the case n — 0. 



The boundary conditions (126) and (127) translate into 

or or 



o„ iujir^nt = rbn + ico^h^ne = rc„ = 



at r = re- 



(158) 
(159) 



The viscous dissipation rate per unit volume is 

2AtS2 + ^b(V • ^t)^ (160) 
and the time- averaged tidal dissipation rate in the convective region can be evaluated as 

(161) 



= 7r / Dn r^ dr, 



where 



dan A„ 



-|- n{n + 1) 



CLn , dbn 

h r—— 

r dr 
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dr 3 

+ {n - l)n{n +l){n + 2) (|6„|2 + |c„|2) 



dCfi 
dr 



(162) 



5.3. Method of solution 

We truncate the system of equations at even order L, so that o„ = 6„ = c„ = (i„ = 
for n > L, and thereby obtain a large but finite number of linear, inhomogeneous ordinary 
differential equations constituting a two-point boundary-value problem. The order of the 
system is in fact 3L. 



5.3.1. Principal method 

Our principal method of solution makes use of a Chebyshev pseudospectral approach 
(e.g., Boyd 2001) as used by, e.g., Rieutord et al. (2001) in their study of free inertial waves 
in an incompressible fiuid contained in a spherical annulus. This method is well suited to 
the problem as it provides spectral accuracy, superior to any finite-difference method with a 
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similar number of radial nodes, and naturally supplies additional resolution in the boundary 
layer just outside the sohd core. 

We introduce the Chebyshev coordinate x such that —l<x<l and 

r=fl^Vo+fii^Vb, (163) 



and reorganize the differential equations by defining a new set of dependent variables 

Xi = C2i-i, Yi = a2i, Zi = 1 < i < ^, (164) 

using equation (157) to substitute for wherever it appears. For each value of i we then 
obtain an equation involving the quantities 

d'X, dX, dY, dY,_, 

dx^ dx dx dx 

from equation (156). Differentiating equation (157) with respect to r, we obtain an equation 
involving the quantities 

f • ^" ^'''^ 

Finally, eliminating dn between equations (153) and (155), we obtain an equation involving 
the quantities 

dXi d^Yi dYi d?Z, dZ, dXi+, 

dx dx'^ dx dx'^ dx dx 

The functions Xj, 1^, and Zj, for the appropriate values of i, are represented by their values 
at the Gauss-Lobatto collocation nodes 

Xj^cos(^-^^, 0<j<iV, (168) 

where N is the Chebyshev truncation order. Each differential equation is then represented 
at each of the internal collocation nodes, using the first and second Chebyshev collocation 
derivative matrices. At the boundary nodes, the boundary conditions are represented in- 
stead, again using the Chebyshev collocation derivative where needed. 

The result of this procedure is a block-tridiagonal system of linear, inhomogeneous 
algebraic equations. Each block is of dimension 3(X+1) and relates to the functions c„_i, a„, 
and hn for a single value of n. There are L blocks along the diagonal, and blocks immediately 
above and below the diagonal, which originate from the Coriolis terms in the equation of 
motion. The blocks are dense owing to the structure of the Chebyshev collocation derivative 
matrix. We solve the system by a standard procedure for block-tridiagonal matrices, which 
involves 0{LN^) operations and has a memory requirement 0{LN'^), although this can be 
reduced to 0{N'^) if ample temporary storage is available on disk. 



-41 - 



5.3.2. Alternative method 

Our alternative method is based on a direct integration of the differential equations 
using a fifth-order Runge-Kutta method with adaptive stepsize control. This method does 
not work when the viscosity is very small, because the equations are then very stiff, but 
we use it to verify the results obtained by the principal method in an accessible parameter 
regime. 

We recast the system as a set of 3L first-order differential equations for the functions 
a2i, 02i, C2i-i, d2i, vr —, ur , 1 < ^ < — , (169) 

where v — /x/p is the kinematic viscosity. To eliminate dan/dr and d'^an/dr'^ from the 
equations, we use equation (157) to express these quantities in terms of 6„ and dbn/dr. 
To start the integration at r = rb we must guess 3L/2 initial values in addition to those 
specified by the boundary conditions, and there are 3L/2 corresponding boundary conditions 
to satisfy at r = Vc- As the problem is linear, a conventional shooting method is not required. 
Instead, we integrate the equations once from Tc to rb using null initial values, to obtain a 
'particular solution'. We then integrate the equations 3L/2 times, omitting the forcing terms, 
setting a different one of the initial values to unity each time (the others being null) . This 
approach generates a set of 'complementary functions' that spans the space of admissible 
initial values. The desired solution is then the particular solution plus a linear combination 
of the complementary functions; the coefficients of the complementary functions are found 
by imposing the boundary conditions at rb, requiring the inversion of a square matrix of 
dimension 3L/2. 



5.4. Computation of the Hough functions 

In order to determine the amplitudes with which the Hough waves are excited in the 
radiative region, according to equation (124), the numerical solution at the boundary r = rb 
must be projected on to the Hough functions. In the case of shellular rotation, the operator 
C given by equation (60) reduces to Laplace's tidal operator, the eigenfunctions of which 
have often been expanded in associated Legendre polynomials (e.g., Longuet-Higgins 1968 
and references therein) . We make use of the identity 

P-(cos 0) = cul^C [^n^r(cos 0) + 2(cos 0) + ^(cos 0)] , (170) 
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valid for n > 2, where 
1 



n{n + 1) + fm — f 



(n — l)n + fm 



+ 



n+l 



> + l)(n + 2) + fm 



(171) 

Bn = - /r:TTf ^ (172) 

n{n + 1) + jm 

and / = 2^2(1) Equation (170) is readily verified by direct substitution. It becomes 
indeterminate in the case m ^ when / = —n{n + l)/m for integers n > m, owing to 
the existence of a null eigenfunction of C corresponding to a toroidal mode, as discussed in 
Section 3.6. 

To solve the eigenvalue problem Cw = Xw we consider the equivalent problem {u}'^^^C)^^w 
Aw, where A = (uj'^^-^X)^^ . In cases when C has a null eigenvalue, we consider the restriction 
of this equation to the subspacc of functions orthogonal to the null eigenfunction. Expanding 
w in normalized associated Legendre polynomials. 



w 



J]e„P-(cose), 



where the sum is over even integers n > 2, we obtain from equation (170) 

A-nGn + Bn-ien-2 + -Sn+ie„+2 — Ae„. 



(173) 
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We truncate at the same order n — N as for the solution in the convective region, and 
thereby obtain an eigenvalue problem involving a tridiagonal, real symmetric matrix. 
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It is clear from this representation that, in the case of a non-rotating planet (/ = 0), the 



matrix is diagonal and the eigenvalues are the diagonal elements [n{n + 1)]^ , in agreement 
with equation (69). More generally, the eigensolutions of equation (175) are readily obtained 
by standard numerical methods. The inner product in equation (124) is easily computed in 
the Legendre polynomial representation, because 
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(cos^) 



^6„,P,"^(cos^) 



sin 9 dO = o^^n- 



(176) 
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To determine the particular solution in the transition region, rather than working with 
equation (119), we project equations (112)-(114) for the particular solution on to spherical 
harmonics in the manner of equation (147), leading to (in an obvious notation) 



-n(n + l)za;6(f) -'^(2^+^ 



+20 



n 



l)(n + l)g„ci^2i + n{n + 2)qn+ic^li = -n{n + 1) 



ip) 



(p) 



(177) 
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dlnr^ 

(n - l)(n + l)g„6?2i + n{n + 2)qr,+ib^l, - 2im^lc^^ = 0, 



l-f(r^a„)-n(n + l)6i^) = 0, 

ar 



(178) 
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where a„ = —iujC,rSne- It is then straightforward to solve algebraically for the coefficients 
dn^ defining the desired function W^P\e). When £ = 2, the only non-zero coefficients are 
0,2, b^2 \ cf^ (when m < 2), 4^-*, d^2 \ ^i^^- course of obtaining this solution, a 

singularity arises if the coefficient of Cn^ in equation (178) vanishes. This gives rise to a 
'toroidal mode resonance' as discussed in Section 4.4. The only such resonances for ^ = 2 
tidal forcing in a planet with shellular rotation occur at cu — — Q/3 (when m = 2), or at 
cu — — Q and a) = — (when m — 1). 



6. NUMERICAL RESULTS 

6.1. Model pairameters 

We present results for a uniformly rotating planet with a solid core of fractional ra- 
dius Tc/Ri = 0.2 and a convective-radiative boundary of fractional radius r\jRi = 0.9. 
These values are suggested by the models of short-period extrasolar planets constructed by 
Bodenheimer, Lin, & Mardling (2001). 

The kinematic shear viscosity u is assumed to be uniform and is parametrized by the 
dimensionless Ekman number 

while the bulk viscosity, which we found to be unimportant, is set to zero. Current models 
(Guillot et al. 2003) suggest that the microscopic viscosity in the interior of Jupiter is of 
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the order of 10~^cm^s~^, corresponding to Ek ~ 10~^^. The eddy viscosity associated with 
turbulent convection is much larger, corresponding to Ek ~ 10~^. However, the characteris- 
tic timescale of the eddies is a year or longer, making the convection inefficient at damping 
tidal disturbances with typical periods of a few days. Therefore the estimate Ek ~ 10"''' 
should probably be reduced by several orders of magnitude, perhaps using the prescription of 
Goldreich & Keeley (1977). In this initial study, we make no attempt to model the detailed 
variations of the effective viscosity, but investigate the systematic behavior of the solutions 
as the uniform kinematic viscosity is reduced into the regime Ek <^ 1. 

Almost all the results reported here arc obtained by the principal method, using a stan- 
dard numerical resolution oi L = N = 100, for which each calculation takes approximately 
30 seconds on a modern workstation. As discussed below, this resolution allows numerically 
converged results to be obtained at Ekman numbers at least as small as 10"^. We have also 
carried out a small number of calculations for isolated forcing frequencies at a very high 
resolution of L = N = 500, each taking several hours on a workstation, but allowing access 
to yet smaller Ekman numbers. 

We have computed the tidal dissipation rate by viscous dissipation in the convective 
region and via the emission of Hough waves at the convective-radiative boundary for forcing 
frequencies in the range —3Q < Co < Sri, comfortably spanning the spectrum of incrtial 
waves. Only the most important case of m = 2 forcing was considered. To investigate the 
effect of a small but non-zero viscosity, the calculations were carried out for Ek = 10~^, 
10-^ 10-^ and 10"^ 

In the case of the uniformly rotating planet, the calculations were also repeated omitting 
certain terms in the governing equations, in order to test the validity of earlier and more 
tractable approaches to the problem. In the traditional approximation (e.g.. Chapman & 
Lindzen 1970) the latitudinal component of the angular velocity is neglected in computing 
the Coriolis force. This amounts to neglecting the QsinO terms in equations (138) and (140), 
while the Hough modes are unaffected. The traditional approximation, where valid, is useful 
because it allows the linearized equations to be solved by separation of variables, although 
we do not make use of that property here. In the no-Coriolis approximation the Coriolis 
force is completely neglected and the Hough modes reduce to spherical harmonics. This is 
equivalent to considering a non-rotating planet except that the Doppler shift of the forcing 
frequency is taken into account. 
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6.2. Dimensionless dissipation rates 

It is convenient to express the tidal dissipation rate in a dimensionless form. The natural 
unit for viscous dissipation is 

Use of this unit means that the dimensionless viscous dissipation rate D^isc/Uvisc depends 
only on cD/Q and Ek (in principle also on rc/ Ri and Vh/Ri if these are varied). We note that 
£)visc captures the viscous dissipation of the equilibrium tide in addition to the damping of 
any inertial waves that are excited as a dynamical tide, and Dyisc is therefore non-zero even 
under the no-Coriolis approximation. 

To calculate the dissipation rate via Hough waves, we note that the waves carry both 
energy and angular momentum fluxes, their ratio being the angular pattern speed uj/m. 
When the waves damp, their angular momentum is deposited locally in the planet. Part of 
the energy of the waves is then required in order to change the spin of the planet, while the 
remainder is dissipated locally. A simple calculation shows that the dissipation rate is 

^Hough = (^^) F, (182) 

where F is the energy flux in equation (124) (factors of e and ordering subscripts being 
omitted) . 

The dissipation rate via the emission of Hough waves is naturally expressed in terms of 
the unit 

t^Hough = /HoughtA'isc, (183) 

where 



is a model-dependent, dimensionless quantity. The first factor in /nough, being the ratio of the 
density at the convective-radiative boundary to the central density of the polytropic model, is 
highly uncertain, being possibly in the range 10~^ — 10~^ for short-period extrasolar planets, 
according to the models of spherically irradiated and tidally heated planets by Bodenhcimer 
et al. (2001). The second factor involves the quantity "D, which is also very uncertain but 
is raised only to a weak power. For planets that are distant from their host stars, such 
as Jupiter, /nough is probably so small that tidal dissipation via Hough waves is entirely 
negligible. 



According to the conventional theory of the equihbrium tide, in which the efficiency of 
tidal dissipation is parametrized by a quality factor Q (e.g., Goldreich & Soter 1966), the 
tidal dissipation rate in a body of negligible rigidity would be 



GMlR\ 
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where /c2 is the second-order Love number, and 
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is a dimensionless number of order e^. We can therefore express the numerically determined 
tidal dissipation rates in terms of effective Q- values given by 



For Jupiter, using the volumetric radius for Ri and the value ~ 0.38 from Gavrilov & 
Zharkov (1977), we obtain /g ^ 0.18; when Jupiter is spun down to a period of three days, 
illustrative of a short-period extrasolar planet, we obtain /g ^ 0.0033. 

We note that Goldreich & Soter (1966) assume that k2 = 3/2, as for a homogeneous 
body, and therefore the effective Q-values appropriate for use with their formulae are ap- 
proximately four times larger than the values we quote below. 



Fig. 2 shows the viscous dissipation rate as a function of forcing frequency, for the 
standard model. The four panels correspond to Ekman numbers 10"*^, 10~^, 10~^, and 
10^^. The quantity plotted is the base-10 logarithm of {D^isc/ Uvisc)\^/ i^l, which is directly 
related to the effective Q- value. The dotted line in panel (d) corresponds to a fiducial value of 
Qvisc — 10^ for Jupiter, or approximately 50 times higher for a short-period extrasolar planet. 
Fig. 3 shows the dissipation rate via Hough waves for the same model. Although the values 
obtained are much larger, they must be scaled down by the small and uncertain factor /nough 
before they can be compared with Fig. 2. Figs 4-7 show the equivalent results obtained when 
either the no-Coriolis approximation or the traditional approximation is adopted. 

It is clear from Fig. 2 that the results become more complicated, and the computations 
more challenging, as the viscosity is reduced. The numerical convergence of the solutions is 
examined in two resolution studies in Figs 8 and 9. The first of these (Fig. 8) is carried out at 
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6.3. Standard model and veiriations 
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a forcing frequency cu — Q, within the spectrum of inertial waves but not associated with an 
obvious resonant feature. The figure shows that higher resolution is required at lower Ekman 
numbers, and that the dissipation rate can be either underestimated or overestimated if the 
resolution is insufficient. The standard resolution L = N = 100 clearly gives numerically 
converged results for all the Ekman numbers considered in Fig. 2. The second resolution 
study (Fig. 9) is carried out at a forcing frequency a) = 0.489^2, which lies near the peak 
of the largest resonant feature seen in Fig. 2. In this case the standard resolution is just 
adequate to capture the peak dissipation rate at Ek — 10"^ with a small error. 

Fig. 10 shows, as a dotted line, an expanded view of part of panel (d) of Fig. 2, the 
graph of viscous dissipation rate for Ek = 10~^. The solid line, for comparison, conies from 
a calculation at Ek — 10~^, for which the increased numerical resolution L — N — 200 was 
necessary. 

Finally, Fig. 11 illustrates the spatial structure of the velocity field induced by tidal 
forcing in the convective region. This calculation, at Ek — 10~^, is performed at the very 
high resolution ol L — N — 500, and the forcing frequency Cj — —1.1181 Q, is chosen to be 
close to the peak of a large resonance. 

6.4. The Coriolis effect and the traditional approximation 

We begin by comparing the results for the uniformly rotating planet in the full model 
(Figs 2 and 3) with those obtained under the no-Coriolis approximation (Figs 4 and 5). 
It is clear that, when the Coriolis force is included, a number of resonant features are su- 
perimposed on the otherwise smooth variation of the viscous dissipation rate with forcing 
frequency. As expected, the resonances are restricted to the interval — 2 < a)/Q < 2 cor- 
responding to the spectrum of inertial waves in the convective region. For Ek = 10"^ the 
inertial waves show up as two resonant peaks in the dissipation rate. However, as Ek is 
reduced further, a host of resonances come into play. This trend is consistent with the idea 
that the spectrum of inertial waves is dense or continuous in the absence of viscosity, while 
it consists of a discrete set of damped modes in the presence of viscosity. When Ek = 10^^ 
the discrete nature of the modes is clearly evident, but when Ek — 10~^ the individual 
resonances can barely be discerned, except for the most prominent examples, and a resonant 
response is almost guaranteed for any forcing frequency within the spectrum of inertial waves. 
Under the no-Coriolis approximation the convective regions of the planet have no wavelike 
response and no dynamical tide is excited. The slow, large-scale motion associated with the 
equilibrium tide is damped by viscosity, leading to a dissipation rate that vanishes linearly 
with the viscosity. The same dependence does not occur in the full model for frequencies 
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within the spectrum of inertial waves. This behaviour reflects the fact that the characteristic 
spatial scale of the tidally forced inertial waves diminishes as the viscosity is reduced. Rather 
as occurs in hydrodynamic turbulence, the dissipation rate then has a non-trivial finite limit 
as the viscosity tends to zero. 

The dissipation rate via Hough waves is differently normalized, owing to the model- 
dependent prcfactor /nough- If this factor is of the order of 10~^ — 10~^, as may occur in 
short-period cxtrasolar planets, the Hough dissipation rate could dominate over the viscous 
dissipation rate. A comparison of Figs 2 and 3 indicates that the Hough dissipation rate 
responds to the incrtial-modc resonances when the full Coriolis effect is used. This result 
reflects the fact that the Hough waves are excited partly by the pressure of the inertial waves 
at the convective-radiative boundary. These resonant features become stronger and sharper 
as the Ekman number is reduced, but otherwise the Hough dissipation rate is essentially 
independent of the viscosity. The feature at Cj — — Q/3 is the toroidal mode resonance 
mentioned in Section 5.4. In the neighborhood of this forcing frequency, a global toroidal 
mode is excited in the radiative region, contrary to the assumptions of this paper, and the 
Hough dissipation rate cannot be trusted. The no-Coriolis approximation fails to reproduce 
any of these features but, generally, significantly overestimates the Hough dissipation rate. 
This discrepancy may be attributable to the fact that, depending on their frequency, the 
Hough modes are confined by the Coriolis force to propagate only in certain ranges of latitude, 
unlike the g modes in a non-rotating planet. 

Now comparing Figs 2 and 3 with Figs 6 and 7, we find that, within the inertial region, 
the traditional approximation gives a totally inaccurate representation of the tidal dissipa- 
tion rate. While the traditional approximation does allow for a resonant response in the 
convective region, the details are quite wrong and the viscous dissipation rate is generally 
overestimated by about two orders of magnitude. This difference may occur because of a 
artificial decoupling between the angular and radial structures of the response that occurs 
when the traditional approximation is employed. This decoupling implies that, when the 
conditions for resonance are met, a global amplification of the response occurs. We conclude 
that, within the spectrum of inertial waves, the traditional approximation is, if anything, 
even worse than the no-Coriolis approximation, and should therefore never be used for iner- 
tial waves in gaseous giant planets. 

6.5. The low- viscosity limit 

Returning to the full model, we see clearly from the expanded view in Fig. 10 that 
the viscous dissipation rate does not simply vanish linearly with the viscosity. Instead, as 
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the Ekman number is reduced (in this case by a factor of 10) the resonant peaks become 
more numerous, sharper, and taller. The frequency-averaged viscous dissipation rate for 
the interval shown in Fig. 10 is very nearly independent of the viscosity, decreasing by 
approximately four per cent when the Ekman number is reduced by a factor of 10. In fact, 
the frequency sampling in Fig. 10 is such that the resonant peaks for the case Ek = 10~^ 
are not perfectly captured, and this inaccuracy may be sufficient to account for the four per 
cent reduction. 

This finding suggests what may be a very important result: although the effective Q- 
valuc derived from the dissipation of inertial waves at very low Ekman number may be an 
exceedingly complicated function of the forcing frequency, a robust average Q-value may be 
obtained that is asymptotically independent of the viscosity as Ek — > 0. In Appendix A 
we present a toy model that may explain this effect. The toy model, based on a simplified 
Cartesian geometry, has the advantage of possessing an analytical solution from which it can 
be shown that the frequency- averaged dissipation rate is exactly independent of the viscosity. 
We suggest that this property may be typical of situations in which forcing is applied to 
systems that exhibit a dense or continuous spectrum in the absence of dissipation. 

It is conceivable that the dissipation rate, averaged locally in frequency, may be es- 
sentially independent of the details of the dissipation mechanism, whether it is due to a 
Navier-Stokes viscosity, damping by turbulent convection, or nonlinear wave damping. We 
are far from having demonstrated such a result. Nevertheless, our results suggest that a 
Q-value of the order of 10^ can be obtained for Jupiter or Saturn through this mechanism. 
The value for a synchronized short-period extrasolar planet would be some 50 times greater 
owing to its slower rotation. 

Fig. 11 illustrates the spatial structure of the tidal response in the convective region. 
In line with work by Rieutord et al. (2001) on free inertial waves in an incompressible fluid 
contained in a spherical annulus, we flnd that the disturbance is localized near rays, which 
are the characteristics of the spatially hyperbolic equations governing inertial waves, and are 
seen to reflect many times from the inner and outer boundaries. In a uniformly rotating 
planet the rays are straight and a simple relation exists between the wave frequency and 
the angle made by the rays with the vertical. A mild near-singularity occurs when the rays 
approach the axis of rotation. The overall wave pattern does not appear as neat as that 
in the calculations by Rieutord et al. (2001), perhaps because the density is non-uniform, 
but more probably because of an interference effect between the equilibrium and dynamical 
tides. 
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6.6. Effects of differential rotation 

We have also made a very preliminary investigation of the possible effects of differential 
rotation. Unfortunately, without a reliable theory of how differential rotation is generated 
by convection and other processes in giant planets and stars, the choice of angular velocity 
profile is unconstrained. We considered the case of a linear profile Q{r) of angular velocity in 
the convective region. There are clear similarities between the graphs of viscous dissipation 
rate versus forcing frequency (which we omit here) and those shown in Fig. 2 for a uniformly 
rotating planet. Even with a differential rotation of 20 per cent, however, the resonant 
response is substantially richer and stronger. The tidal response also looks qualitatively 
similar to that shown in Fig. 11, although the rays are curved in a differentially rotating 
planet. 

We emphasize that we have not considered the possible effects of corotation resonances. 
While such resonances in the radiative region may absorb outgoing Hough waves, as noted by 
Goldreich & Nicholson (1989), there may also be a direct tidal torque and energy dissipation 
associated with tidal forcing at the corotation resonance, as occurs in differentially rotating 
disks (e.g., Goldreich & Tremaine 1980). There is also the possibility of an angular velocity 
that depends on latitude. The numerical method used in this paper is not well suited to 
such a situation, unless Q depends in a sufficiently simple way on 9 that the couplings 
between different spherical-harmonic components of the velocity field are strongly restricted. 
Therefore there remains much to be investigated with regard to tides in differentially rotating 
planets and stars. 

7. COMPARISON WITH PREVIOUS WORK 

The previous work most closely related to the present paper is that by Savonije & 
Papaloizou (1997) and Papaloizou & Savonije (1997). The former paper presents a numerical 
analysis of the linearized response of a uniformly rotating high-mass star to i = m = 2 tidal 
forcing. Unlike the planet considered in our calculations, the star is predominantly radiative 
and has only a small convective core. In comparison with our study, Savonije & Papaloizou 
(1997) take a more direct numerical approach to the full linearized equations, including 
non-adiabatic effects and a small viscosity, but retaining only terms of first order in the 
angular velocity. Perhaps because the convective core is small, they find that main effects 
of rotation occur through the resonant excitation of Hough modes in the radiative region, 
rather than inertial modes in the convective region. Indeed, they find a number of toroidal 
mode resonances for negative forcing frequencies, whereas precisely one such resonance is 
allowed in our analysis. Our calculation is more dedicated to the exacting task of resolving 
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the tidally forced inertial waves in the convective region and examining more carefully the 
limit of small viscosity. The reduction of the basic equations that we carry out allows us 
to solve the resulting equations much more efficiently, and the numerical method that we 
employ is particularly well suited to this problem. 

The accompanying paper by Papaloizou & Savonije (1997) provides an asymptotic anal- 
ysis of low-frequency tidal forcing in the same type of star. There is much in common between 
their concisely expressed analysis and the present paper. We have attempted to give a fuller 
exposition of the problem, formally justifying the approximations involved and allowing for 
the important possibility of differential rotation. We have thereby generalized the derivation 
of reduced equations for inertial modes and Hough modes. Furthermore, our analysis should 
be useful for future direct numerical simulations of tides, as it clarifies the role of the anelas- 
tic and Cowling approximations. These approximations are satisfied by the dynamical tide 
but not by the equilibrium tide (in Papaloizou & Savonije 1997 the effect of self-gravitation 
on the equilibrium tide is neglected). 

Related work for non-rotating stars and planets has also been published by Goldreich 
& Nicholson (1989), Lubow et al. (1997), and Goodman & Dickson (1998), who arrive at 
asymptotic expressions for the tidal dissipation rate through the emission of g modes at a 
convective-radiativc interface. There is a close formal relation between these expressions 
and our equation (124). However, the expression can only be evaluated, in the case of a 
rotating star or planet, if the inertial wave problem in the convective region has been solved 
numerically. In this connection it would be interesting to rc-cxaminc the problem of tidal 
forcing in rotating solar-type stars with convective envelopes (cf. Terquem et al. 1998; 
Goodman & Dickson 1998) taking into account either uniform or differential rotation. In 
the latter case the angular velocity profile of the solar convective envelope, inferred from 
helioseismic studies, could be used as a guide. 

8. SUMMARY AND CONCLUSIONS 

In this paper we have revisited the classical problem of the response of a giant planet to 
low-frequency tidal forcing. Much of our general analysis of this problem applies equally to 
tidal forcing in stars. The novel feature of our approach is that we take into account the slow 
and possibly non-uniform rotation of the planet. (Here 'slow' rotation means that we include 
the full Coriolis force but not the centrifugal distortion of the planet.) Convective regions of 
the planet then support low-frequency inertial waves with an intricate spatial structure and 
a rich frequency spectrum, while radiative regions support (generalized) Hough modes, some 
of which may be regarded as g modes of high radial order, modified by the (differential) 
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rotation. We have argued that, in many cases of interest, such as that of an extrasolar 
planet that orbits close to its host star, the effective tidal forcing frequency lies within the 
spectrum of inertial waves in the convcctivc regions of the planet, and a resonant response 
can be expected. As a result, the rate of tidal dissipation may be greatly enhanced relative 
to a model that neglects the rotation of the planet. The dissipation occurs both through the 
viscous or turbulent dissipation of the inertial waves in the convective region and through the 
emission of generalized Hough waves that propagate through the radiative envelope towards 
the surface, where they presumably damp. Enhancement of the tidal dissipation rate implies 
a more rapid synchronization of the planet's spin with its orbit, a faster circularization of 
the orbit, and a more intense heating of the planet, which may lead in turn to inflation and 
even Roche- lobe overflow (Gu et al. 2003). 

We have presented a systematic asymptotic analysis of the linearized response of a 
differentially rotating planet to tidal forcing with a frequency that is small compared to the 
dynamical frequency of the planet. The response separates naturally into an equilibrium 
tide, which represents a large-scale, quasi-hydrostatic distortion of the planet in the imposed 
tidal potential, and a dynamical tide, which constitutes a mostly wavelike correction. We 
obtain the reduced system of equations governing the dynamical tide in convective and 
radiative regions separately, and explain the asymptotic matching procedure between the 
two solutions. 

In convective regions, which we model as being adiabatically stratified, the dynamical 
tide takes the form of an indirectly forced inertial wave confined in a spherical annulus. The 
reduced equations are intrinsically two-dimensional and require a numerical solution. More- 
over, for effective forcing frequencies within the dense or continuous spectrum of free inertial 
waves, the problem is mathematically ill-posed unless viscosity is included. A question of 
particular interest is how the dissipation rate varies with the viscosity in the limit that the 
viscosity tends to zero. 

In radiative regions, the dynamical tide involves predominantly horizontal motions and 
takes the form of a wave (or evanescent disturbance) with a short radial wavelength. A sepa- 
ration of variables is possible, and we flnd that the angular structure of the wave is governed 
by a generalization of Laplace's tidal operator. We call the resulting solutions generalized 
Hough modes and discuss some mathematical properties of the operator concerned. 

Our analysis clarifies the role in the tidal problem of certain well known approximations 
to the equations of fluid dynamics. We flnd that the so-called traditional approximation, 
in which the latitudinal component of the angular velocity is neglected, is not applicable to 
inertial waves in convective regions and should never be used for this purpose. In fact, it 
gives highly misleading results for the tidal dissipation rate, perhaps worse than neglecting 
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the Coriolis force altogether. On the positive side, we argue that the dynamical tide satisfies 
both the anelastic and Cowling approximations, which eliminate acoustic waves and self- 
gravitation, respectively, from the problem. However, the equilibrium tide satisfies neither 
of these approximations. These considerations are important not only for analytical or semi- 
analytical studies but also for direct numerical simulations of tidal forcing. Such simulations 
will be needed in order to determine properly how the dynamical tide interacts with turbulent 
convection, as well as the role of nonlinearity in the tidal problem. It is highly convenient 
and appropriate to simulate convection with a numerical method based on the anelastic 
and Cowling approximations. Our analysis indicates that the dynamical tide can also be 
captured within this approach, and shows how the subtle, indirect forcing of the dynamical 
tide by the equilibrium tide can be achieved through the inertial terms in the equation of 
motion. 

We have presented the results of full numerical calculations of the tidal response for 
an idealized model of a giant planet that is predominantly convective but also contains 
a solid core and has a thin radiative envelope. The numerical method is suited to the 
case of 'shellular' rotation in which the angular velocity is independent of latitude, and we 
have focused mainly on the case of uniform rotation. High-resolution calculations, using 
a pseudospectral method, are required to access the physically interesting regime of small 
Ekman number (small viscosity) and reveal the intricate spatial structure of the inertial 
waves while properly resolving the dissipative structures. We have calculated the tidal 
dissipation rate, both by viscous dissipation in the convective region and via the emission of 
Hough waves, as a function of the tidal forcing frequency, for the important case of a tidal 
potential proportional to the £ — m = 2 solid harmonic, and for forcing frequencies that 
span the spectrum of inertial waves. Wc find that the viscous dissipation rate is strongly 
enhanced, relative to a calculation in which the Coriolis force is neglected, in a number of 
inertial-mode resonant peaks that become more numerous, sharper, and taller as the Ekman 
number is reduced. As a result, the viscous dissipation rate does not vanish linearly with the 
viscosity. Depending on the physical conditions at the convective-radiative boundary, the 
Hough dissipation rate may possibly exceed the viscous dissipation rate, and also responds 
to the inertial-mode resonances. A single 'toroidal mode resonance' is also possible, in which 
residual forcing excites a large-scale mode (related to a Rossby wave), rather than a short- 
wavelength response, in the radiative region. When the planet rotates differentially, the 
inertial-mode resonances are yet stronger and more numerous, and the viscous dissipation 
rate is further enhanced. 

Examination of the spatial structure of the tidal response in the convective region shows, 
in line with work by Rieutord et al. (2001) on free inertial waves in an incompressible fluid 
contained in a spherical annulus, that the disturbance is localized near rays, which are the 
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characteristics of the spatially hyperbolic equations governing inertial waves, and are seen 
to reflect many times from the inner and outer boundaries. The rays are straight or curved 
depending on whether the planet rotates uniformly or differentially. 

The dissipation rate associated with inertial waves at very low Ekman number is, in 
principle, a highly erratic function of the forcing frequency (Fig. 2). As we will discuss in 
future work, when realistic scenarios for tidal evolution are considered, there are various 
reasons why it may be more appropriate to apply a smoothed version of this 'fractal' curve 
rather than taking every resonant peak literally. We have presented both numerical evidence 
(Fig. 10), and an analytical demonstration based on a toy model (Appendix A), that the 
frequency-averaged tidal dissipation rate associated with inertial waves may be asymptoti- 
cally independent of the viscosity in the limit of small Ekman number. If correct, this result 
is most important because the Ekman number based on the microscopic viscosity in the 
interior of Jupiter is exceptionally small and beyond the range of any numerical calculation. 
Even the eddy viscosity based on turbulent convection gives rise to a very small Ekman num- 
ber, which should probably be further reduced owing to the mismatch between the typical 
tidal forcing frequency and the characteristic timescale of the convective motion. 

The resulting tidal dissipation rates are not adequately represented by a constant Q- 
value as is commonly adopted in parametrized models. One reason for this is that the 
dissipation rate associated with inertial waves scales naturally with the spin frequency of the 
planet in a way that is not captured in the constant-Q model. Nevertheless, the effective 
Q-values obtained are of the order of 10^ for Jupiter or Saturn, or some 50 times greater for 
a synchronized short-period extrasolar planet owing to its slower rotation. 

These values are probably adequate to explain the historical evolution and current state 
of the Galilean satellites (e.g., Peale 1999). We consider the alternative model of loannou 
& Lindzen (1993b) improbable as it requires the interior of Jupiter to be subadiabatically 
stratified, in contradiction of current models (Guillot et al. 2003) which suggest that it has 
a minuscule superadiabatic gradient owing to the high efficiency with which convection can 
transport the required heat flux. 

In the case of short-period extrasolar planets, a Q-value of the order of 10^ following 
synchronization is probably sufficient to explain the circularization of their orbits (cf. Marcy 
et al. 1997). The additional route of dissipation via the emission of Hough waves in an 
outer radiative layer may play a role here, although this mechanism is subject to significant 
uncertainties. The dissipation of inertial waves in the convective region provides a deep 
source of heating that can inffate the planet, as may have been observed in the case of the 
transiting planet HD 209458 b (Brown et al. 2001). 
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In conclusion, we have shown that it is both important and feasible, although numeri- 
cally challenging, to include the effects of rotation when studying the response of a planet or 
star to tidal forcing. Our general analysis lays the ground work for future numerical studies 
including more realistic interior models and, possibly, the effects of differential rotation or 
nonlinearity. In our preliminary numerical investigation we have demonstrated that a robust 
enhancement of tidal dissipation results from the inclusion of rotational effects, which can 
account for the rapid evolution of tidally interacting systems. There remains much of interest 
to be explored in this problem. 
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A. DISSIPATION OF INERTIAL WAVES IN A TOY MODEL 

In this appendix we present a toy model that explains some aspects of the viscous dissi- 
pation of forced inertial waves. Although the model does not capture many of the subtleties 
of the tidal-forcing problem, it has the advantage of possessing an analytical solution. 

We consider an incompressible fluid of uniform density p and kinematic viscosity v. The 
fluid initially rotates with uniform angular velocity fi, with a pressure gradient balancing 
the centrifugal force and any gravitational force. 

A body force / per unit mass excites small disturbances satisfying the linearized equa- 
tions (written in the rotating frame) 

^ + 212 X w = — Vp + i/V^w + /, (Al) 
at p 

V • u = 0, (A2) 

where u and p are the Eulerian perturbations of velocity and pressure. 

The fluid is contained in a square box, < x < L and < 2; < L, referred to Cartesian 
coordinates {x,y,z) such that 17 = Oe^. For the purposes of illustration, we assume that 
all quantities are independent of and adopt the boundary conditions 

Ux^Uy^ ^ = (A3) 
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on X — 0, L and 

t^ = ^ = «. = (A4) 

oz oz 

on ^ = 0, L. Apart from the first condition on Uy, these can be interpreted as meaning that 
the boundaries are rigid and stress-free. 

We adopt units of mass, length, and time such that p = 1, L = tt and 2Vt = 1. For a 
periodic forcing, all quantities may be assumed to have the form 



= Re [u{x, z) e"^'^*] , Uy^Re [v{x, z) e"^'^*] , 
= Re [w{x, z) e-^'^*] , p = Re [^l^{x, z) e"''^*] , 



U = Re [fix, z) e-'^'] , = Re [h{x, z) e"^'^*] . (A5) 

We assume for simplicity that fy — 0. 

The linearized equations then reduce to 

-iuu-v = -d^i) + v{dl + dl)u + f, (A6) 

-iujv + u = u{dl + dl)v, (A7) 

-iojw = -d^ip + u{dl + dl)w + h, (A8) 

d^u + d,w = 0, (A9) 



subject to 
and 



u — V — dxW — at X — 0,7r (AlO) 



dzU — dzV — w — Q at 2; = 0, TT. (All) 
For free modes of oscillation (f — h — 0), solutions exist of the form 

u = Umn sin COS nz, V = siu mx cos nz, 

W = Wjnn COS TUX siu HZ, 1p = Iprnn COS mX COS HZ, (Al2) 

with m and n non-negative integers, not both equal to zero. The dispersion relation of these 
viscously damped inertial modes is 

= (A13) 

where 

uj = uj + ii'{m'^ + n^). (A14) 
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In the absence of viscosity, the eigenfrequencies are real and dense in the interval — 1 < a; < 1. 
Viscous damping causes the eigenfrequencies to move below the real axis, increasingly so for 
modes of higher 'order'. 

We return to the forced problem where lo is real and prescribed. The force components 
/ and h can be expanded in Fourier series (the sums are over non-negative integers m, n) 

/ = ^ fmn sin mx COS nz, 

h — hmn COS mx sin nz. (A15) 

The solution to the forced problem is then 

u — sin mx cos nz, 

V — Vmn sin mx cos nz, 

w = Wmn COS mx sm nz, 

■0 = 'ipmn COS mx COS nz, (A16) 



with 



iujn{nfmn - mh 



{m? + n'^)Cj'^ — n 



2 ' 



(A18) 



_ n(l -a;2)/i^„-ma;Vmn /.^^x 

Of particular interest is the average rate of viscous dissipation per unit volume, 

D=— / {-u-iyV'^u)dxdz. (A21) 

TT^ Jo Jo 

This evaluates to 



and can then be related to the forcing in the form 

8 |(m^ + n^)a;^ — n 
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The variation of D with u and v is very comphcated. As — > for fixed u, D eventually 
tends linearly to zero, unless ui is one of the eigenvalues of the inviscid problem (which 
form a dense set of measure zero), in which case D tends to infinity. However, an 'average' 
dissipation rate associated with the spectrum of inertial waves can be identified by integrating 
D over all positive (or negative) frequencies, with the simple result 



which is exactly independent of the viscosity. 

Fig. 12 illustrates how D varies with uj for several values of z/, for a simple forcing of 
the form fmn = ^/mn'^, hmn = 0. (In fact, for generic forcing, = Oil/min?) for large 
m, n.) As seen in the tidal-forcing problem, the dissipation is enhanced in the neighborhood 
of inertial-mode resonances in the interval —l<uj<l. As u is reduced the resonant peaks 
become more numerous, sharper, and taller, but the integrated dissipation rate remains 
exactly the same. In fact, the contribution of each mode to the integrated dissipation rate is 
independent of the viscosity, and the strongest resonances, which are with modes of lowest 
'order', retain their identity and importance as v is reduced. 




(A24) 
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Fig. 1. — Tidal forcing frequencies and the spectrum of inertial waves in a uniformly rotating 
planet. The ratio of the (spin) angular velocity of the planet to the mean motion of the orbit 
is plotted against the ratio of the effective forcing frequency to the angular velocity of the 
planet, for components 2-5 of the tidal potential. The dotted hnes indicate the extent of the 
spectrum of inertial waves. A short-period extrasolar planet being spun down towards the 
synchronized state descends the tracks from near the top of the diagram to the dashed line. 
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Fig. 2. — Tidal dissipation rate by viscosity in a uniformly rotating planet, according to the 
full model. The dissipation rate is scaled for comparison with equation (187). Panels (a), (b), 
(c), and (d) are for Ekman numbers 10~^, 10~^, 10~^, and 10"^, respectively. The spectrum 
of inertial waves corresponds to the interval [—2,2]. The dotted line in panel (d) indicates 
the dissipation rate in Jupiter corresponding to an effective quality factor Qvisc = 10^. 
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Fig. 3. — Tidal dissipation rate via Hough waves in a uniformly rotating planet, according 
to the full model. The dissipation rate is scaled for comparison with equation (187). Panels 
(a), (b), (c), and (d) are for Ekman numbers 10~^, 10~^, 10~^, and 10~^, respectively. 
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Fig. 4. — Tidal dissipation rate by viscosity in a uniformly rotating planet, according to the 
no-Coriolis approximation. The parameters are as in Fig. 2. 
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Fig. 5. — Tidal dissipation rate via Hough waves in a uniformly rotating planet, according 
to the no-Coriolis approximation. The parameters are as in Fig. 3. 
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Fig. 6. — Tidal dissipation rate by viscosity in a uniformly rotating planet, according to the 
traditional approximation. The parameters are as in Fig. 2. 
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Fig. 7. — Tidal dissipation rate via Hough waves in a uniformly rotating planet, according 
to the traditional approximation. The parameters are as in Fig. 3. 
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Fig. 8. — Resolution study away from resonance. The viscous dissipation rate for a forcing 
frequency cu — Q,, not associated with an obvious resonant feature, is plotted versus the 
Ekman number for numerical resolutions L — N — 500 (solid line), 200 (dotted hne), 100 
(dashed line), 50 (dash-dotted line), 20 (dash-triple-dotted line), and 10 (long-dashed line). 
The three highest resolutions give results that are almost indistinguishable by eye. The 
alternative numerical method reproduces the curve for Ek > 10~^ to very high accuracy, 
but cannot be used for smaller Ek. 
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Fig. 9. — Resolution study on resonance. As for Fig. 8, but for a forcing frequency a) = 
0.489 associated with the largest resonant feature seen in Fig. 2. 
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Fig. 10. — Expanded view of the graph of viscous dissipation rate in the full model, for 
forcing frequencies close to uj = Q. The dotted line corresponds to Ek = 10^^ and is 
equivalent to panel (d) of Fig. 2. The solid line corresponds to Ek = 10~® and is calculated 
at a higher resolution oi L — N — 200. The peaks of the curve are not perfectly captured at 
the frequency samphng of Aoj/Q — 0.002. Although the viscosities differ by a factor of ten, 
the frequency-averaged dissipation rates over this interval differ by at most a few per cent. 



Fig. 11. — High-resolution calculation {L = N = 500, Ek = 10~^) of the tidal response of a 
uniformly rotating planet. The root-mean-square velocity of the total tide (equilibrium plus 
dynamical) is plotted in a meridional slice through the convective region only. The velocity 
scale is linear, black representing zero. The forcing frequency is chosen to be near the peak 
of an inertial-mode resonance. 
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Fig. 12. — Variation of the viscous dissipation rate with forcing frequency, in the toy model 
for forced inertial waves (Appendix A). Panels (a), (b), (c), and (d) are for dimensionless 
viscosities v — 10~^, 10~^, 10~^, and 10~^, respectively. The spectrum of inertial waves 
corresponds to — 1 < a; < 1, and the graphs are symmetric about uj — Q. 



